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ABSTRACT 


The  diffractive  optic  imaging  spectrometer,  DOIS,  is  a  high  resolution,  compact, 
economical,  rugged,  programmable,  multi-spectral  imager.  The  design  implements  a 
conventional  CCD  camera  and  emerging  diffractive  optical  element  (DOE)  technology  in  an 
elegant  configuration,  adding  spectroscopy  capabilities  to  current  imaging  systems  [Lyons 
1995], 

One  limitation  of  DOEs,  also  known  as  zone  plate  lenses,  is  abundant  chromatic 
aberration.  DOIS  exploits  this  typically  unwanted  effect,  utilizing  a  DOE  to  perform  the 
imaging  and  provide  the  dispersion  necessary  to  separate  a  multi-spectral  target  into 
separate  spectral  images.  The  CCD  is  stepped  or  scanned  along  the  optical  axis  recording  a 
series  of  these  spectral  images.  This  process  is  referred  to  as  diffractive  spectral  sectioning. 

Under  this  dissertation,  three-dimensional  spectral/spatial  DOE  imaging  theory  was 
developed  to  describe  and  predict  the  system’s  performance.  The  theory  was  implemented 
in  a  software  model  to  simulate  DOIS  image  cubes.  A  visible  spectrum  DOIS  prototype 
was  designed,  fabricated  and  characterized.  The  system's  incoherent  point  spread  function 
was  theoretically  modeled  and  experimentally  determined.  To  verify  the  simulations,  the 
prototype's  performance  was  demonstrated  with  a  variety  of  known  targets  and  compared 
to  simulated  image  cubes.  To  reconstruct  the  three-dimensional  object  cubes,  various 
deconvolution  algorithms;  nearest  neighbor,  inverse  filtering  and  constrained  iterative 
deconvolution,  were  developed  and  applied  to  both  computer  generated  and  experimentally 
measured  image  cubes.  The  best  results  were  obtained  using  an  SVD  inverse  Fourier 
deconvolution  algorithm  with  regularization  for  noise  suppression.  The  results  demonstrate 

a  resolving  power  greater  than  288  (X/AX=577nm/2nm).  Finally,  three  additional  DOIS 
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designs  are  presented  as  suggestions  for  future  work,  including  a  configuration  with  no 
moving  parts  which  records  the  entire  3D  image  cube  in  one  "snapshot". 

DOIS  is  a  practical  image  spectrometer  that  can  be  built  to  operate  at  ultraviolet, 
visible  or  infrared  wavelengths  for  applications  in  surveillance,  remote  sensing,  medical 
imaging,  law  enforcement,  environmental  monitoring,  and  laser  counter  intelligence. 
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1.  INTRODUCTION 


The  field  of  multi-  and  hyper-spectral  imaging,  also  known  as  Imaging 
Spectrometry,  has  a  long  history  [Goetz  1995]  and  has  been  receiving  increased  attention 
for  several  years.  As  discussed  elsewhere  [Descour  1995],  Imaging  Spectrometry  adds  the 
ability  to  examine  the  spectral  distribution  of  two-dimensional  scenes  to  the  fundamental 
power  of  imaging  systems.  The  availability  of  known  spectral  radiance,  reflectance  and 
absorption  curves  coupled  with  an  imaging  spectrometer  allows  identification  and 
classification  of  targets  with  an  accuracy  and  resolution  previously  unknown. 

An  entire  chapter  of  the  IR  Handbook  [Wolfe  1989]  is  dedicated  to  the  spectral 
properties  of  natural  sources.  Recorded  spectra  can  be  compared  to  these  and  other  known 
spectra  to  identify  the  material  composition  of  a  target.  A  spectrometer  is  an  instrument  for 
resolving  and  recording  these  spectra.  The  combination  of  an  imager  and  a  spectrometer,  an 
imaging  spectrometer,  provides  a  conventional  spatial  image  as  well  as  the  spectral  content 
of  each  pixel,  forming  a  three-dimensional  spatial/spectral  image  cube.  The  result  is  a 
complete  system  capable  of  target  detection,  classification  and  identification. 

Image  spectrometers  can  be  built  to  operate  in  the  ultra-violet  (UV),  visible  or 
infrared  (IR)  spectrum  and  measure  spectral  reflectance,  emission  and/or  absorption. 
Spectra  can  be  used  as  a  "fingerprint"  to  determine  exactly  what  material  or  mix  of  materials 
is  present.  Whether  the  data  is  from  the  paint  of  a  military  tank,  the  plume  of  a  missile,  the 
ground  in  the  desert  or  the  air  we  breathe,  spectral  information  can  lead  to  recognition  of  a 
detected  target  or  even  detect  the  presence  of  unwanted  materials. 

1.1  Image  Spectrometry 

The  goal  of  an  image  spectrometer  is  to  obtain  a  set  of  spectral  images  of  a  target  to 
form  an  image  cube.  Each  image  represents  a  spectral  band  or  channel,  defined  by  a  central 


wavelength,  Xc,  and  its  spectral  bandwidth,  AX.  The  total  number  of  channels  is  equal  to 
the  entire  wavelength  range  divided  by  the  channel  bandwidth. 

#  of  channels  =  -  Xmin)/AX  (1.1) 

Common  image  spectrometers  include  three  components:  a  dispersion  element,  an 
imager  and  a  detection  device  [Wolfe  1989].  Examples  include  a  lens  and  color  filter  wheel 
in  front  of  a  two-dimensional  detector  array,  or  a  single  pixel  aperture  illuminating  a 
diffraction  grating,  dispersing  the  light  onto  a  linear  detector  array.  The  detector  records 
one  pixel  at  all  wavelengths  which  is  then  scanned  in  the  x  and  y  directions  to  form  an 
image. 

1.2  Figures  of  Merit 

When  designing  or  comparing  imaging  spectrometers  it  is  customary  to  discuss  the 
spectral  resolution  of  the  system  [Wolfe  1994],  Spectral  lines  usually  have  some 
predetermined  shape,  such  as  a  Gaussian  or  Lorentzian.  The  line  width,  sometimes  defined 
as  the  spectral  slit  width,  is  usually  specified  as  the  full  width  at  half  maximum  (FWHM), 

as  shown  in  Figure  1.1,  noted  as  AX. 

spectral  resolution  =  A> l  (1-2) 

The  relative  resolution  is  the  spectral  resolution  divided  by  the  center  wavelength. 

relative  resolution  =  AX/XC  (1.3) 

Since  bigger  is  always  better,  the  resolving  power  is  defined  as  the  reciprocal  of  the  relative 
resolution.  resolving  power  =  Xj 'AX  (1.4) 
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Figure  1.1:  Typical  spectral  line  used  to  define  spectral  resolution. 

1.3  Common  Image  Spectrometers 

The  most  simplistic  image  spectrometer  is  a  standard  camera  (lens  and  detector) 
with  various  line-pass  color  filters.  Typically  the  filters  are  placed  on  a  wheel  that  is  rotated 
so  that  they  are  in  the  optical  path  one  at  a  time.  Figure  1 .2  shows  that  the  data  recorded  is  a 
series  of  two-dimensional  (x,y)  images  stepped  in  wavelength.  The  number  of  filters 
determines  how  completely  the  spectra  is  sampled  and  the  spectral  resolution  of  each  filter 
is  the  spectral  resolution  of  the  spectrometer.  However,  the  system's  size,  complexity  and 
cost  are  directly  related  to  the  number  and  quality  of  the  filters,  and  the  materials  needed  for 
some  filters  have  limited  availability.  Additionally,  once  designed  and  built  it  is  inflexible; 
interest  in  additional  wavelengths  will  require  an  upgrade. 


Figure  1.2:  Three-dimensional  image  cube  for  a  color  filter  wheel  image 
spectrometer. 
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Another  popular  design  is  the  grating  spectrometer.  A  narrow,  rectangular  entrance 
slit  is  imaged  onto  a  linear  diffraction  grating  and  then  on  to  a  two-dimensional  detector. 
The  grating's  dispersion  is  perpendicular  to  the  long  axis  of  the  slit.  The  detector  records  a 
spatial  image  in  one  direction  and  the  spectra  along  the  other.  The  data  is  a  set  of 
spatial/spectral  images,  (x,A)  or  (y,A),  that  are  recorded  as  the  aperture  is  scanned  along  the 
scene  in  a  pushbroom  manner  to  form  an  image  shown  in  the  image  cube  of  Figure  1.3.  Its 
advantage  is  that  there  is  a  continuous  spectra  from  the  diffraction  grating,  so  the  resolution 
is  determined  by  the  optical  magnification  and  detector  size  (parameters  that  a  designer 
should  have  full  control  over).  However,  in  order  to  view  a  conventional  (x,y)  image  for 
target  acquisition  and  tracking,  it  requires  a  scanning  mechanism  which  can  be  sensitive  to 
vibration  and  motion  of  the  target. 


y 


Figure  1.3:  Grating  image  spectrometer  three-dimensional  image  cube. 
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A  third  design  is  a  Michelson  interferometer  setup  with  an  oscillating  mirror  in  one 
path.  As  the  mirror  moves,  the  interference  term  varies  as  the  cosine  of  the  oscillation  with 

a  frequency  of  Infk  for  a  monochromatic  source.  A  source  of  a  different  wavelength  or 

frequency  will  provide  an  interference  pattern  with  a  different  frequency,  2kI\2.  If  the  two 

are  used  at  the  same  time  then  the  pattern  will  be  the  sum  of  two  cosines,  each  modified  by 
its  amplitude.  With  a  multi-spectral  source,  the  interferogram  obtained  as  the  mirror  moves 
is  the  sum  of  a  collection  of  monochromatic  interference  patterns,  each  with  its  own 
amplitude.  The  interference  pattern  is  recorded  over  time,  correlated  to  the  optical  path 
difference  caused  by  the  mirror,  and  Fourier  transformed  to  yield  the  spectrum.  A  CCD 
focal  plane  array  can  be  placed  at  the  detection  plane  to  simultaneously  record  the  frequency 
spectrum  of  multiple  pixels,  forming  an  image  once  the  transform  is  performed,  as  shown 
in  Figure  1 .4.  This  technique  provides  very  high  resolution  data  however  it  is  sensitive  to 
vibrations  and  target  motion  during  recording. 


Figure  1.4:  Three-dimensional  data  cubes  for  a  Fourier  Transform  image 
spectrometer. 

1.4  Spectrotomography  Techniques 

Spectrotomography  is  a  new  branch  of  optical  tomography  where  the  three- 
dimensional  data  cube  is  integrated  along  various  axes  forming  two-dimensional 
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projections  of  the  3D  data.  These  projections,  containing  integral  data  on  the  spectral  and 
spatial  properties  of  the  object  of  interest,  are  recorded  then  processed  with  computer 
tomography  algorithms  to  reconstruct  the  3D  data  cube.  Researchers  use  several  techniques 
with  various  dispersion,  scanning,  and  reconstruction  methods  for  applying  computer 
tomography  to  imaging  spectrometry.  A  major  advantage  of  these  systems  is  high 
throughput  efficiency.  The  projections  integrate  all  possible  photons  and  therefore  no  light 
is  discarded,  which  is  useful  in  low  light  scenarios.  Although,  the  recorded  data  requires 
much  processing  to  form  a  recognizable  image  cube.  Also,  much  emphasis  has  been  placed 
on  the  future  ability  to  build  a  system  with  no  moving  parts  and  with  a  single  snap-shot 
approach  for  moving  targets  [Descour  1995], 

1.5  DOIS 

The  diffractive  optic  image  spectrometer  design  is  a  spectro-tomographic  technique 
which  integrates  over  a  narrowband  of  the  image  cube.  Shown  in  Figure  1.5,  the  DOE 
provides  both  the  spatial  imaging  and  the  spectral  dispersion  in  this  image  spectrometer 
[Hinnrichs  1995,  Lyons  1995], 

A  panchromatic  CCD  camera  is  used  for  detection.  This  camera,  mounted  parallel  to 
the  DOE,  is  stepped  along  the  optical  axis,  z .  Each  detector  location  z  corresponds  to  a 
specific  spectral  channel  of  the  target’s  image.  As  an  example,  Figure  1.6  lists  the 
calculated  image  distance  z{  and  magnifications  for  particular  spectral  slices  of  an  object 

located  2 fd  from  the  DOE,  where  fd  is  the  focal  length  at  the  DOE  design  wavelength  kd. 

The  camera  is  interfaced  to  a  frame  grabber  board  within  a  computer  for  image  capture, 
analysis,  display  and  post  detection  processing.  At  each  step  along  the  optical  axis,  both  an 
image  frame  and  the  corresponding  detector  location  is  recorded. 
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Figure  1.5:  Three-dimensional  Spectral  imaging  with  a  diffractive  optical 
element  (DOE). 
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Figure  1.6:  The  spectral  image  distances  zt  and  magnifications  for  an  object 
at  z=2fd. 

The  images  captured  by  the  CCD  are  a  superposition  of  one  channel’s  infocus 
image  and  defocused  images  of  the  surrounding  channels.  These  defocused  components 
cause  blurring  of  the  desired  image. 

Like  a  color  filter,  DOIS  integrates  over  a  bandpass,  dictated  by  the  DOE’s 
incoherent  point  spread  function.  However,  DOIS  can  provide  greater  spectral  resolution 
by  scanning  along  the  z  axis,  shifting  the  center  of  the  bandpass  and  sampling  the  spectra  at 
several  spectral  points  for  every  point  that  the  filter  could  sample.  Once  sampled,  digital 
image  restoration  techniques  can  be  applied  to  remove  the  blurred  information  and 
reconstruct  the  original  spectra. 
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The  DOE  dispersion  is  so  great  that  each  image  is  an  integral  of  data  in  only  a 
narrow  spectral  band,  not  the  entire  spectrum.  This  leads  to  a  system  that  is  useful  both 
with  and  without  post-detection  processing.  Without  processing  DOIS  provides  spectral 
images  with  spectral  resolution  that  rivals  color  filter  wheels.  The  spectral  resolution  is 
defined  by  the  region  of  spectral  integration  measured  to  be  as  good  as  1.5nm  for  spatially 
unresolved  targets  [Chapter  7],  Higher  resolving  power  is  achievable  with  DOIS  by 
applying  digital  image  restoration  algorithms  which  improve  the  resolution  by  eliminating 
spectral  and  spatial  blur  [Chapter  8]. 

1.6  Dissertation’s  Objectives 

Parallel  to  the  efforts  under  this  dissertation,  Michele  Hinnrichs  from  Pacific 
Advanced  Technology  (PAT)  conceived  and  patented  the  mechanical  design  of  the 
diffractive  optic  image  spectrometer  described  in  this  dissertation  [Hinnrichs  1995]. 
Although  an  IR  system  had  been  built  and  demonstrated  on  real-world  targets,  there  was  a 
lack  of  analysis  theory  and  understanding  of  spectral  imaging  with  a  diffractive  optic.  This 
theory  was  needed  for  analytical  modeling  and  simulation  of  the  image  spectrometer 
design,  and  required  to  develop  object  reconstruction  algorithms  for  imaging  resolved 
targets  with  high  spectral  and  spatial  resolution. 

The  objectives  of  this  dissertation  were  to  determine  the  theory  behind  diffractive 
spectral  imaging,  implement  the  theory  in  a  software  model  to  simulate  DOIS  image  cubes, 
experimentally  verify  the  simulated  results  with  a  hardware  prototype  imaging  known 
spectral/spatial  targets,  and  develop  and  apply  object  reconstruction  algorithms  to  the 
simulated  and  experimentally  measured  image  cubes.  Additionally,  an  improved  design 
was  conceived  to  provide  diffractive  optic  image  spectrometry  with  constant  magnification 
[Lyons  1997]. 
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In  the  next  chapter,  the  capabilities  of  the  diffractive  optic  image  spectrometer  are 
first  demonstrated.  A  visible  spectral  range  DOIS  prototype  imaged  an  AppleColor™  high 
resolution  RGB  computer  monitor  with  the  words  AFMC  ROME  LAB  spelled  out  in  blue, 
green  and  red  phosphor  respectively.  Although  a  much  more  analytical  characterization  is 
presented  later  in  Chapter  7,  this  experiment  is  useful  for  discussion  and  introduction  to 
DOIS's  operation. 

Since  the  main  component  is  a  diffractive  optical  element  (DOE)  providing  both 
imaging  and  dispersion  for  this  image  spectrometer,  Chapter  3  details  some  important 
theory,  design  and  characteristic  equations,  fabrication  techniques  and  modeling  schemes 
of  the  Fresnel  Zone  Plate  which  is  the  DOE  used  in  DOIS. 

The  theory  of  image  formation  with  the  diffractive  optic  is  developed  in  Chapter  4. 
The  three-dimensional  image  cube  is  collected  by  recording  a  set  of  images  taken  at 
sequential  focal  planes.  By  analogy  with  physical  sectioning  techniques,  this  method  is 
called  diffractive  spectral  sectioning.  For  convenience,  Figure  1.7  lists  the  six  Cartesian 


Coordinate  systems  used  throughout  this  dissertation. 


object  space  =>  obj(x0,y0,h0)  at  z, 
Gaussian  image  space  =>  o(x,y,z) 
image  space  =>  i(x,,y,,z,) 
isotomic  image  space  =>  i  (x/  ,y-  ,z-) 
isotomic  Gaussian  image  space  =>  o'  (x'  ,y'  ,z') 

reconstructed  object  space  =>  objr{x0 , y0 ,  X0)  at  z„ 


Figure  1.7:  The  six  Cartesian  Coordinate  systems  used  in  DOIS. 

The  DOIS  prototype  was  fully  characterized,  first  as  individual  components  and 
then  as  a  complete  system.  This  characterization  is  presented  in  Chapter  5. 

Chapter  6  presents  both  computer  generated  (CG)  and  experimentally  determined 
image  sets  of  known  targets.  As  you  will  see,  the  linear  system  theory  of  Chapter  4 
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accurately  predicts  output  images  formed  by  DOIS.  The  CG  images  are  very  similar  to 
those  determined  experimentally. 

A  multi-spectral,  multi-spatial  four  target  testbed  was  designed  and  built  to  further 
test  the  performance  of  DOIS.  The  four  targets  include  a  Tungsten-Halogen  greybody 
source,  a  multi-line  Mercury  source  and  two  Helium  Neon  lasers,  a  632.8nm  HeNe  and  a 
542nm  GreNe.  Chapter  7  summarizes  the  testbed's  known  spectral  characteristics  and 
includes  a  collection  of  spectral  images  and  spectral  radiance  curves  measured  with  DOIS. 
Appendix  A  is  a  catalog  of  almost  600  images  of  the  four  targets  that  were  recorded  with 
the  DOIS  prototype.  The  most  compelling  results  are  the  images  and  spectral  plots  of  the 
Mercury  doublet  at  577nm  and  579nm.  Without  any  post-detection  processing  this  doublet 
is  resolved,  demonstrating  a  spectral  resolution  of  less  than  2nm! 

Chapter  8  presents  the  object  reconstruction  algorithms  to  generate  in-focus  spectral 
images,  spectral  distributions,  and  full  three-dimensional  representations.  Since  not  all 
applications  require  the  same  spectral  resolution,  processing  techniques  which  fall  in  three 
general  categories  are  presented,  each  providing  different  amounts  of  deblurring  at  various 
computational  expense;  nearest  neighbor,  inverse  filtering  and  constrained  iterative 
deconvolution. 

The  algorithms  were  implemented  in  Mathematica  and  applied  to  reconstruct  both 
CG  and  experimentally  measured  objects.  The  Nearest  Neighbor  algorithm  estimates  any 
spectral  image  with  only  three  "snapshots"  and  minimal  computational  expense.  However, 
greater  resolution  can  be  obtained  using  an  SVD  inverse  Fourier  deconvolution  with 
regularization  for  noise  suppression.  For  completeness,  constrained  iterative  deconvolution 
is  presented,  although  in  practice  this  method  has  a  relaxation  function  which  must  be 
optimized  to  a  particular  application. 

All  of  the  deconvolution  algorithms  provide  improved  resolution  and  remove  the 
effects  of  large  targets,  higher  diffraction  orders  and  stray  light.  The  choice  between 


24 


reconstruction  approaches  will  ultimately  depend  on  the  application  dictating  such  issues  as 
the  required  spectral/spatial  resolution,  data  storage,  computation  time  and  memory. 

Chapter  9  summarizes  the  conclusions  drawn  during  the  course  of  this  research, 
recommending  the  most  effective  approach  to  transition  this  work  into  a  fieldable  image 
spectrometer.  The  advantages  of  DOIS  over  conventional  image  spectrometers  are 
discussed.  Finally,  three  designs  for  future  DOIS  systems  are  presented;  a  dual-waveband 
design  for  simultaneous  spectral  imaging  in  two  spectral  bands,  such  as  the  mid  IR  (3  to 
5pm)  and  far  IR  (8  to  12pm),  and  two  configurations  that  perform  image  spectrometry 
with  no  moving  parts  including  a  design  which  records  the  entire  3D  image  cube 
simultaneously  in  one  "snapshot". 
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2.  DOIS  -  DIFFRACTIVE  OPTIC  IMAGE  SPECTROMETER 


The  principle  shortcoming  of  DOEs  for  polychromatic  imaging  is  abundant 
chromatic  aberration.  This  chapter  presents  a  patented  design  that  exploits  this  typically 
unwanted  effect  to  create  an  image  spectrometer  [Hinnrichs  1995,  Lyons  1995],  Under  this 
dissertation  a  visible  regime  Diffractive  Optic  Image  Spectrometer  (DOIS),  shown  in  Figure 
2.1,  was  designed,  simulated,  fabricated,  demonstrated  and  characterized.  A  DOIS  can  be 
designed  to  operate  at  ultraviolet,  visible  or  infrared  wavelengths  for  multispectral  and 
hyperspectral  imaging  in  medicine,  forensics,  industrial  and  environmental  monitoring,  as 
well  as  military  applications. 

From  Chapter  1,  the  goal  of  an  image  spectrometer  is  to  obtain  a  set  of  spectral 
images  of  a  target,  forming  an  image  cube.  Each  spectral  image  represents  a  spectral  band 

or  channel,  defined  by  a  central  wavelength,  Xc,  and  its  spectral  bandwidth,  AX.  Referring 

to  equation  (1.1),  the  total  number  of  channels  is  equal  to  the  total  wavelength  range 
divided  by  the  channel  bandwidth. 

#  of  channels  =  (Xmax- X^J/AX 

2.1  Basic  Scenario 

A  color  computer  monitor  serves  as  the  polychromatic  source  object;  it  will  be 
referred  to  as  the  target.  The  DOE  provides  both  the  spatial  imaging  and  the  spectral 
dispersion  in  this  image  spectrometer.  A  conventional  monochromatic  CCD  camera  is  used 
to  record  the  image  cube.  This  camera,  mounted  parallel  to  the  DOE,  is  stepped  along  the 
optical  axis,  z .  Each  detector  location,  z,  corresponds  to  a  specific  spectral  channel  of  the 
image  cube.  At  each  step  both  an  image  frame  and  the  corresponding  detector  location  is 
recorded. 
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Figure  2.1:  DOIS  -  Diffractive  Optic  Image  Spectrometer  basic  scenario. 


The  camera  is  interfaced  to  a  frame  grabber  board  within  a  computer  for  image 
capture,  analysis,  display  and  post  detection  processing.  The  images  captured  by  the  CCD 
are  a  superposition  of  one  channel’s  infocus  image  and  defocused  images  of  the 
surrounding  channels.  These  defocused  components  cause  blurring  of  the  desired  image. 
Post-detection,  digital  image  restoration  algorithms  [Gonzalez  1987]  are  used  to  remove  the 
unwanted  blurred  components,  improving  the  imager’s  spatial  and  spectral  resolution 
[Chapter  8]. 

The  diffractive  optic  1st  order  imaging  equation  relating  spectra  and  the  object  and 
image  distances  ,  s0  and  is  given  in  equation  (2.1)  [see  Chapter  3  for  the  derivation]: 

J__  J _ 1_  _  __A _ 1_ 

sol  fx  So  Kfd  So  (2.1) 

where  fx  is  the  focal  length  at  wavelength  X,  calculated  from  fd  and  Xd  the  DOE  design  focal 
length  and  wavelength. 
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solving  for  wavelength: 

if  « then 


f  i  i  N 
^  =  ^dfd - ^  — 

^  _  ^dfd  _  K dfd 
sa  fx 


(2.2) 

(2.3) 


Figure  2.2  shows  the  use  of  equation  (2.2)  for  the  scenario  in  Figure  2.1. 

With  the  target  relatively  close  to  DOIS  as  in  the  laboratory  demonstration,  it  is  best 
to  keep  the  DOE  fixed  and  the  detector  stepping,  since  an  accurate  target  to  DOE  distance  s0 
is  required.  The  practical  field  implementation  would  be  designed  for  far  range  targets, 
s0~°o.  At  infinite  conjugates  the  object  distance  s0  isn’t  sensitive  to  changes  in  the  lens 
location,  the  DOE  can  be  scanned  and  the  detector  can  be  kept  stationary,  a  preferable 
scenario. 


s  =  7  ft 


L  =  20  cm 


=  588  nm 


sn  =  20.63 cm 

K  -  Kfl — +— 


sn  =  25.3  cm 

K  =  Kfi— +- 


sn  =  29.73cm 


^3  “  ^dfd  - + 


1  1 


A,  =  588  *  20 
A[  =  625 nm 


1  1 
20.63  213.36 


A,  =588*2of— 3— +  — - 

^  l  er  o  n  -i  *1  ■ 


25.3  213.36 


Aj  =  520 nm 


A3  =  588*2of — —  +  — - — 
V  29.75  213.36 

A,  =  450 nm 


Figure  2.2:  Algorithm  for  calculating  wavelength  from  image  distance,  s-. 

When  a  target  such  as  AFMC  ROME  LAB  spelled  out  on  a  black  background  in 
blue,  green  and  red  respectively,  is  imaged  with  DOIS,  each  color  or  word  will  come  to 

focus  in  a  different  position  zx.  The  chromatic  focal  shift,  Af  =  %  is  so  severe  that 


the  separation  between  a  blue  AFMC  and  a  red  LAB  is  almost  10  cm.  This  extreme  fall  off 
with  color  provides  this  image  spectrometer  with  a  FWHM  spectral  bandwidth  ranging 
from  1  to  1.5  nm  throughout  the  visible  spectrum  and  over  100  possible  spectral  channels 
without  post  detection  image  processing.  Once  a  series  of  images  is  captured,  the  image 
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cube  can  be  displayed  as  series  of  spectral  images,  a  normalized  spectral  radiance  curve 
(intensity  verses  wavelength)  for  each  pixel  of  interest  or  a  two-dimensional  spatial/spectral 

cross-section,  xzx  or  yzx,  of  the  image  cube.  The  stepping  can  be  controlled  to  tune  to  one 
particular  channel,  to  step  through  a  narrowband,  wideband  or  the  entire  visible  range  in 
equal  increments  of  dz,  d?i,  or  any  variation  of  user  controlled  steps.  With  proper  control,  a 

continuous  video  can  be  recorded  matching  the  frame  and  scan  rates,  to  calibrate  the 
wavelength  to  time  or  frame  number. 

2.2  Experimental  Demonstration  with  a  RGB  monitor 

The  Image  Spectrometer  of  Figure  2.1  was  built  and  demonstrated  at  Rome 
Laboratory’s  Photonic  Center  and  is  described  below.  This  monitor  demonstration  is 
presented  to  facilitate  an  understanding  of  the  basic  operation  and  performance  of  DOIS.  A 
much  more  analytical  testing  of  the  characteristics,  performance  and  limitations  of  this 
Image  Spectrometer  is  found  in  the  following  chapters;  Chapter  5.  Prototype 
Characterization  and  Chapter  7.  Prototype  demonstration  with  four  targets. 

The  DOE  is  a  5  cm  diameter,  BK-7,  2  level  phase  grating  with  a  design  focal 

length,  fa  =  20  cm  at  the  design  wavelength,  Aj  =  588  nm.  The  CCD  camera  is  a  SONY 

XC-75,  monochrome  video  camera.  The  sensing  area  is  1/2  inch  with  768  (H)  x  494  (V) 

* 

pixel  elements.  Each  pixel  is  8.4  pm  x  9.8  pm.  The  approximate  pixel  size  of  10  pm  x  10 
pm  will  be  used  for  simplicity  and  to  account  for  dead  space  between  pixels.  The  camera 
was  connected  to  a  framegrabber  board  in  a  Macintosh  PowerPC  8100/100AV.  The  images 
were  recorded  and  displayed  with  MacPhase  by  Otter  Solution,  an  image  capture,  analysis 
and  visualization  software  package. 


29 


The  target  simulation  monitor  used  in  this  demonstration  is  an  AppleColor™  High 
Resolution  RGB  monitor.  Each  pixel  of  an  RGB  monitor  is  a  cluster  of  three  phosphor 
dots:  a  red,  a  green  and  a  blue.  All  colors  displayed  on  the  monitor  are  a  combination  of 
these  three  basic  colors.  For  instance,  white  is  equal  amounts  of  all  three,  bright  yellow  is 
equal  amounts  of  red  and  green,  pink  is  blue  and  red,  and  turquoise  is  blue  and  green. 
Adobe  Photoshop™  is  a  software  program  which  allows  precise  selection  of  text  colors 
and  was  used  to  create  the  display  target.  The  words  AFMC  ROME  LAB  were  spelled  out 
in  Photoshop™  on  a  black  background.  The  color  picker  was  used  to  select  100%  blue  font 
color  for  the  AFMC,  and  100%  green  and  red  for  the  ROME  and  LAB  respectively.  Figure 
2.3  is  an  example  of  the  color  picker  window  in  Adobe  Photoshop™  shown  here  with 
100%  blue  for  the  AFMC  text. 


Ored  A 
0  green 
255  blue 


Figure  2.3:  Example  of  the  color  picker  window  in  Adobe  Photoshop™  for 
the  blue  AFMC  target. 

Figure  2.4  shows  the  relative  spectral  emission  curves  of  the  red,  green  &  blue 
phosphors.  The  general  wideband  spectral  emission  curves  for  the  three  phosphors;  blue 
P22B,  green  P22G  and  red  P22R  [Fink  1989].  Experimental  characterization  of  the 
monitor  s  spectral  radiance  was  performed  with  an  Oriel  grating  monochrometer  and 
Newport  power  meter  .  This  characterization  confirmed  the  wide  band  spectral  emissions 
of  the  blue  and  green  phosphors.  However,  the  red  phosphor  emitted  in  narrow  lines,  with 
two  strong  red  lines,  625  and  617  nm,  a  few  orange  lines  around  594  nm  and  even  a  few 
blips  in  the  green  region.  Since  the  red  emission  differed  so  greatly  from  the  textbook 
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P22R  spectral  curve,  both  are  plotted  in  Figure  2.4.  Unsuccessful  attempts  were  made  to 
verify  the  spectral  radiance  of  the  monitor  with  the  manufacturer.  The  line  nature  of  the  red 
could  indeed  be  from  the  P22R  phosphor  or  a  possible  coating  or  filter  on  the  monitor 
screen,  and  the  lines  in  the  green  region  could  be  from  misalignment  of  the  electron  beam 
accidentally  exciting  the  green  phosphor  dots.  In  any  event,  the  band  and  line  emission  of 
this  monitor  makes  an  excellent  demonstration  target. 

The  horizontal  curve  in  Figure  2.4  is  the  DOIS  minimum  detection  threshold.  DOIS 
is  able  to  detect  a  faint  green,  an  orange  and  the  two  strong  red  lines  of  the  red  phosphor 


and  the  continuous  wideband  emissions  of  the  blue  and  green. 


Monitor  Spectral  Curve 


Figure  2.4:  Phosphor  emission  curves  for  red  (P22R),  green  (P22G),  and 
blue  (P22B)  phosphor,  the  measured  blue  and  green  agree  with  wideband 
textbook  curves  however  the  measured  red  differs  dramatically. 

Spectral  images  of  the  monitor  taken  with  DOIS  are  shown  in  Figure  2.5.  These  are 
raw  (256  x  256)  images,  displayed  without  any  processing.  The  images  depict  the  band 
emission  of  P22B,  P22G  and  the  line  nature  of  P22R.  Scanning  from  400nm  to  650nm  the 
blue  AFMC  is  present  alone  (a)  from  400nm  to  490nm,  at  490  the  green  ROME  cuts  on 
and  both  AFMC  and  ROME  can  be  seen  in  varying  intensities  (b)  until  the  blue  falls  off 
completely  at  540nm.  The  green  ROME  is  continually  present  until  595nm.  The  intensity  of 
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Figure  2.5:  Spectral  images  recorded  with 
DOIS.  The  target  is  text  displayed  on  an 
AppleColor™  high  resolution  RGB  monitor. 

(a)  X  =  452  nm,  blue  AFMC,  mid  blue  focus.  Blurring 
is  from  the  wide  band  emissions  of  P22B. 


(b)  X  -  520  nm,  green  ROME,  mid  green  focus  with  the 
blur  from  the  band  emissions  of  P22G.  Note  that 
corresponding  to  Figure  2.4,  P22B  is  still  emitting. 


(c)  X  =  594  nm,  green  ROME  and  red  LAB.  Both  P22G 
and  P22R  have  slight  emissions  at  this  orange  line. 


(d)  X  =  617  nm,  red  LAB.  First  strong  P22R  line  with  a 
halo  from  625  nm  line. 


(e)  X  =  625,  red  LAB.  Strongest  P22R  line.  Note  that 
P22G  is  completely  cutoff. 
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the  red  LAB  pops  in  and  out  as  the  scan  progresses  from  540nm  to  650nm.  A  very  dim 
LAB  shows  up  at  542nm,  then  again  at  594nm  as  in  (c).  It  is  absent  until  615,  with  a  crisp 
focus  at  617  and  625nm,  finally  disappearing  at  630nm. 


2.2.1  Blue  Background  Demonstration 

The  monitor  demonstration  was  repeated  with  a  blue  background  instead  of  the 
black.  This  simulates  a  scenario  with  a  target  surrounded  by  a  bright  background  (the  sky). 


Figure  2.6:  Images  from  green  ROME  and  red  LAB  on  a  blue  background 
(a)  blue  focus  (b)  green  (c)  red. 

As  shown  above,  DOIS  is  able  to  distinguish  the  targets.  The  defocused  blue  light 
adds  to  the  noise  floor  and  decreases  the  SNR.  Notice  that  the  green  ROME  is  clear  even 
though  the  blue  phosphor  is  still  emitting  in  that  region. 
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2.3  Continuous  Magnification 

Recall  from  Figures  1.5  and  1.6  that  the  magnification  changes  with  wavelength. 
Since  the  CCD  pixel  configuration  is  fixed,  the  image  cube  is  recorded  with  non-uniform 
sampling.  This  creates  problems  in  registering  target  information  between  spectral  slices 
and  limits  the  accuracy  of  image  reconstruction  algorithms. 


2.3.1  Post-detection  image  resampling 

To  compensate  for  the  changing  magnification  the  image  cube  can  be  resampled  to 
approximate  constant  magnification.  Illustrated  in  Figure  2.7,  a  pixel  footprint  or  field-of- 
view  (FOV)  is  defined  at  the  smallest  magnification.  Each  spectral  image  is  then  resampled 
by  integrating  over  the  same  pixel  FOV,  Qpixe]. 


_  \txei  __  integration  area 

llpixel  2  ~  2 

z,  z2 


(2.4) 

(2.5) 


^ 'pixel 


_  lOpmxlOpm  _  integration  area 


20  cm2 


30cm2 


(2.6) 


integration  area  =  225pm2  =  1 ,5pixels 


(2.7) 


In  order  to  form  a  constant  magnification  image  cube  with  each  spectral  image 
having  256x256  spatial  pixels,  the  red  image  at  z=20cm  is  recorded  with  256x256  and  the 
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blue  image  at  z=30cm  will  need  to  be  recorded  with  384x384  pixels,  then  interpolated  to  fit 

the  256x256  spectral  slices  of  the  image  cube. 

Problems  can  arise  with  this  technique.  Resampling  the  detected  image  cube 
requires  averaging  a  non-integral  number  of  pixels.  This  can  cause  blurring  of  edges  and 
inaccurate  pixel  registration.  A  more  accurate  optical  method  of  obtaining  an  image  cube 
with  constant  magnification  is  presented  in  section  2.3.2. 

2.3.2  Optical  designs  for  constant  magnification 

A  second  approach  to  DOIS  imaging  with  constant  magnification  is  a  redesigned 
optical  train.  A  refractive  lens  is  added  between  the  DOE  and  CCD  to  relay  and  magnify  the 
DOE  images  [Lyons  1997].  The  relay  lens  is  the  scanned  element,  and  its  position  dictates 
the  spectrum  of  the  recorded  image. 


constant.  The  DOE  still  performs  the  dispersion  and  forms  an  image  space,  however  a 
refractive,  non-dispersive  lens  is  placed  after  the  DOE  image  space  relaying  the  DOE 
images  to  the  CCD  camera.  The  images  formed  by  this  intermediate  lens  on  the  CCD  will 
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be  erect  and  the  spectra  will  depend  on  the  lens'  position,  /,  and  the  target's  object  distance, 
s0targef  The  following  equations  derive  the  relationship  between  spectra  and  lens  location. 

1  \  ' 

Starting  from  equation  (2.2):  X  =  X dfd  - + - 

^  ^ iDOE  ^target  y 


SiDOE  ~  l  Solens 

(2.8) 

f  y 

o  _  J  lens  Hens 

^  ole  ns  r 

^  Hens  J lens 

(2.9) 

Silens  =  T  —  l 

(2.10) 

.  _  flens(T~l) 

olens  m  j  /> 

i  -  /  -  Jkns 

(2.11) 

T-l-f'ens  T-l-  fUns 

(2.12) 

X  -  X  f  f .  T-l-  flens  1 

A  AdJd  irp  _  j2  _  r\  r  J_  r  rp 

^  J lens  J lens  ^  ^otarget  J 

(2.13) 

The  selected  wavelength  X  is  calculated  in  equation  2.13  from  the  DOE  design  constants 

X/d,  the  tube  length  T,  the  focal  length  of  the  relay  lens  flens,  the  target's  object  distance 
sotarget,  and  the  lens  position  /. 
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zoom  designs. 

The  change  in  magnification  verses  wavelength  for  this  relay  design  is  plotted  in 
Figure  2.9  along  with  that  of  the  original  design.  Notice  that  although  the  magnification 
isn't  constant,  the  change  over  the  entire  spectral  band  is  substantially  reduced  to  less  than 
10%.  To  further  improve  the  changing  magnification,  the  relay  lens  can  be  replaced  with  a 
pair  of  lenses  forming  a  zoom  lens,  shown  in  Figure  2.10.  The  constant  magnification  of 


the  zoom  lens  design  is  shown  as  a  dotted  line  in  Figure  2.9. 


50  "lOO  150  200  250  300  mm 


Figure  2.10:  Zoom  DOIS  design  with  constant  magnification. 
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3.  DIFFRACTIVE  OPTICAL  ELEMENTS 


The  three  main  considerations  in  designing  an  optical  system  are  cost,  weight  and 
performance.  Diffractive  optical  elements  (DOEs)  have  the  potential  to  improve  all  three. 
DOEs  have  a  spectral  dispersion  inversely  proportional  to  wavelength,  whereas  refractive 
optical  elements  have  dispersion  proportional  to  wavelength.  Typically  lens  designers 
combine  these  dispersion  characteristics  to  obtain  an  achromat  [Ward  1971,  Stone  1990], 
However  in  this  system  the  DOE  dispersion  is  used  alone,  and  exploited  to  create  an 
imaging  spectrometer. 

3.1  Fresnel  Zone  Plate 

The  potential  usefulness  of  diffractive  optical  elements  has  been  known  for  years. 
[Miyamoto  1961,  Lesem  1969],  The  problem  with  implementing  them  in  real  systems  has 
been  the  lack  of  a  process  that  can  design  and  reliably  produce  the  elements  to  the  required 
tolerances  necessary  for  high  diffraction  efficiency.  Figure  3.1  (a)  shows  an  example  of  a 
Fresnel  zone  plate  phase  profile  needed  to  achieve  high  efficiency.  The  2n  phase  depth 
corresponds  to  a  material  etch  depth  of  about  2jam  for  midinfrared  radiation  and  lpm  for 
visible.  Several  research  teams  are  exploring  techniques  to  produce  this  continuous  phase 
profile  with  laser  writers  or  greyscale  masks  [Morris  1994,  O'Shea  1994],  However,  it  is 
common  practice  to  approximate  the  continuous  phase  profile  by  quantized  discrete  phase 
levels  [Goodman  1970,  Swanson  1989],  Figures  3.1  (b),  and  (c)  show  the  Fresnel  zone 
plate  phase  profile  quantized  into  two  and  four  phase  levels,  respectively.  The  two  level 
profile  results  in  an  efficiency  of  40%  and  the  four  level  profile  results  in  an  efficiency  of 
81%.  For  most  systems,  it  is  necessary  to  achieve  diffraction  efficiency  of  99%  or  greater, 
this  can  be  achieved  with  16  levels  created  with  4  masks.  Figure  3.2  shows  the  diffraction 
efficiency  as  a  function  of  the  number  of  phase  levels  and  required  masks.  Since  this 
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dissertation  is  only  meant  as  a  proof-of-concept,  40%  is  tolerable  and  a  2  level  phase  plate 
was  fabricated  with  one  mask. 
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Figure  3.2:  Diffraction  efficiency  vs.  number  of  masks  and  phase  levels 
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3.2  DOE  design  and  characteristic  equations 


The  zones  of  the  Fresnel  Zone  Plate  are  defined  when  the  optical  path  difference 

between  adjacent  zones  is  A/2,  causing  constructive  interference,  shown  in  Figure  3.3  and 

characterized  by  equation  3.1  [Hecht  1989],  From  Pythagorean’s  theorem,  the  radius  of  the 
m*  zone,  rm,  is: 


rm  +/ 2  =  f  +  m 1 


(3.1) 


where  m  is  the  zone  counter,  rm  is  the  value  of  the  mth  radius,  A  is  the  design  wavelength, 
and  f  is  the  primary  focal  length,  corresponding  to  the  1st  diffractive  order. 

A  , 

since  m—  «  f 

2  (3.2) 

rm  =  finh 


Solving  for  f: 


(3.3) 
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The  zone  radii  are  calculated  to  focus  a  particular  design  wavelength,  Ad,  at  a  design  focal 


length,  fd. 


rl=fdmXd  (3-4) 

Once  the  DOE  is  fabricated  with  these  fixed  zones  determined  from  the  design  parameters, 
the  focal  position  at  wavelengths  other  than  the  design  wavelength  can  be  found  from 
substituting  equation  (3.4)  into  (3.3): 


f_r2m  =  fdmXd 
mX  mX 


(3.5) 


This  wavelength  dependence  is  seen  as  chromatic  aberration  where  the  effective 
focal  length  is  inversely  proportional  to  wavelength.  Adding  it  to  the  Gaussian  lens  formula 
yields  the  diffractive  optic  1st  order  imaging  equation: 

1  _  1 _ J_  =  ^A _ J_ 

5,.(A)  f(X)  sa  Xdfd  sa  (3.6) 


where  s0  is  the  object  distances  and  s,(A)  is  the  spectral  image  distance. 

A  critical  figure  of  merit  for  this  application  is  the  first  order  diffraction  efficiency, 
Hj.  It  is  dependent  on  the  etch  depth,  the  operational  wavelength  and  the  number  of  phase 

levels  etched  into  the  grating.  The  diffraction  efficiency  can  be  calculated  from  equation 
(3.7)  below,  and  is  plotted  in  Figure  3.4  for  a  2  level  zone  plate  and  a  high  efficiency  16 
level  DOE  [Swanson  1989]. 
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where:  d=Xd/(n-l) 

r|  is  the  diffraction  efficiency 
n  is  the  refractive  index  of  the  substrate 
d  is  the  etch  depth 

m  is  the  diffracted  order  of  interest,  typically  m=l 
Xd  is  the  design  wavelength 
X  is  the  operational  wavelength 

N  is  the  number  of  levels  in  the  grating,  DOIS  prototype  is  2 


The  curve  in  Figure  3.4  is  useful  in  selecting  the  proper  design  wavelength  for 
applications  where  it  is  important  to  have  high  efficiency  through  a  spectral  range  and  not  at 
just  one  wavelength,  such  as  in  DOIS.  Note  that  the  efficiency  falls  off  much  faster  for 

smaller  wavelengths,  therefore  Xd  should  be  selected  near  the  shorter  wavelength  end  of  the 
range  not  simply  the  center  wavelength. 

An  interesting  note  is  that  the  peak  for  the  two  level  grating  is  not  at  X/Xd=l  but  at 

1.3.  This  leads  to  the  conclusion  that  calculation  for  etch  depth,  d,  which  is  formulated  for 
a  continuous  grating,  isn’t  optimal  when  you  have  only  2  levels. 
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3.3  Fabrication  of  Multilevel  Profiles 

Swanson  and  Veldkamp  developed  a  scheme  that  can  accurately  and  reliably 
produce  these  multilevel  diffraction  surfaces  with  diffraction  efficiencies  in  excess  of  90% 
[Veldkamp  1989].  Their  technique  takes  advantage  of  technology  the  integrated  circuit 
industry  has  developed  for  the  miniaturization  of  circuits.  Three  essential  tools  that  were 
developed  are  electron  beam  pattern  generators,  reactive  ion  etchers,  and  mask  aligners.  E- 
beam  pattern  generators  are  capable  of  drawing  binary  amplitude  patterns  with  feature  sizes 
of  0.1  pm  and  positioning  the  features  to  an  even  greater  accuracy.  Reactive  ion  etchers  can 
etch  a  binary  profile  to  depths  of  a  few  microns  with  an  accuracy  on  the  order  of  tens  of 
angstroms.  Mask  aligners  are  used  routinely  to  align  two  patterns  with  an  accuracy  of 
fractions  of  a  micron.  These  are  key  technological  advances  that  make  it  possible  to 
produce  the  high  quality  diffractive  phase  profiles. 

Electron  beam  pattern  generators  produce  masks  that  have  binary  transmittance 
profiles.  A  thin  layer  of  chrome  on  an  optically  flat  quartz  substrate  is  patterned  by  the  e- 
beam  machine.  The  input  to  the  e-beam  pattern  generator  is  a  file  stored  on  a  computer  tape 
and  properly  formatted  for  the  particular  machine.  For  multilevel  diffractive  elements  the 
number  of  phase  levels  in  the  final  diffractive  element  constructed  from  these  masks  is  2n, 
where  n  is  the  number  of  masks.  For  example,  only  four  masks  will  produce  16  phase 
levels,  resulting  in  an  efficiency  of  99%. 

The  binary  amplitude  masks  produced  from  the  pattern  generator  are  actually  low 
efficiency  zone  plates  themselves,  but  are  used  to  construct  the  multilevel  optical  element. 
The  fabrication  process  using  the  first  mask  is  shown  in  Figure  3.5  .  The  optical  surface  on 
which  the  diffractive  profile  is  to  be  etched  is  coated  with  a  layer  of  photoresist.  The  e- 
beam  generated  mask  is  then  placed  over  the  substrate  and  illuminated  with  a  standard  UV 
photoresist  exposure  system.  The  photoresist  is  then  developed,  resulting  in  a  properly 
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patterned  layer  of  photoresist.  The  photoresist  acts  as  an  etchstop  for  the  reactive  ion 
etching. 


Reactive  ion  etching  is  a  process  in  which  an  RF  electric  field  excites  a  gas, 
producing  ions.  The  ions  react  with  the  material  of  the  substrate  and  etch  away  the  surface 
at  a  controlled  rate.  The  reactive  ion  etching  process  is  anisotropic,  so  the  vertical  side  walls 
of  the  discrete  phase  profile  are  retained.  Typical  etch  rates  are  on  the  order  of  100A  to 

o 

200A  per  minute.  As  an  example,  the  required  first-level  etch  depth  for  a  BK-7  substrate  to 
be  used  at  a  wavelength  of  .588  jam  is  1.176|im.  The  necessary  etch  time  is  on  the  order  of 
an  hour.  Numerous  elements  can  be  etched  simultaneously.  After  the  pattern  of  the  first 
mask  has  been  etched  into  the  substrate,  any  residual  photoresist  is  stripped  away. 

For  multilevel  structures  the  same  procedure  is  then  repeated  on  the  optical 
substrate,  only  this  time  the  second  mask  and  etching  to  one-half  the  depth  of  the  first  etch. 
For  the  second  and  subsequent  masks  an  additional  complication  arises.  These  masks  have 
to  be  accurately  aligned  to  the  already  existing  pattern  produced  from  the  first  etch. 
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Commercially  available  mask  aligners  are  capable  of  aligning  two  patterns  to  a  fraction  of  a 
micron.  This  accuracy  is  sufficient  to  retain  diffraction-limited  performance  for  most 
multilevel  structures  designed  to  operate  in  the  visible  and  infrared. 

This  technique  was  used  at  The  Nanofabrication  Facility  at  Cornell  University  to 
fabricate  a  DOE  lenses  for  a  DOIS  prototype  [Lyons  1994]..  The  DOE  is  a  BK-7,  2  level 
Fresnel  Phase  Plate  with  approximately  2200  zones  and  a  minimum  feature  size  of  2.5pm. 


Figure  3.6  is  two  scanning  electron  microscope  (SEM)  photos  of  the  rings. 


3.4  The  Incoherent  Point  Spread  Function 

The  key  parameter  for  a  successful  simulation  and  reconstruction  is  an  accurate  psf, 
h.  Since  the  equivalence  between  DOEs  and  lenses  has  been  established  [Sweatt  1979],  the 
three-dimensional  psf  of  a  lens  can  be  adapted  to  the  DOE.  The  theoretical  diffraction 
limited  and  experimentally  measured  intensity  point  spread  functions  are  presented  below. 


3.4.1  Diffraction  limited  incoherent  psf 

Bom  and  Wolf  describe  the  normalized  three-dimensional  distribution  of  intensity 
in  the  neighborhood  of  focus,  I(v,u),  for  a  rotationally  symmetric  optical  system  as  a 
convergent  series  of  Bessel  functions  with  the  two  equivalent  expressions  [Bom  1989], 
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Figure  3.7  is  a  density  plot  of  equation  (3.8),  in  terms  of  the  dimensionless 
variables  u  and  v,  and  the  scaled  dimensions  r[(im]  and  z[|im]  calculated  from  the  specific 
parameters  of  the  DOIS  prototype  design.  In  three-dimensions  h(v,u)  is  circularly 

symmetric  about  the  optical  axis,  u.  When  evaluated  with  specific  design  parameters,  X,  f, 
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and  aperture  radius,  a,  the  intensity  point  spread  function  will  be  scaled  but  maintain  the 
fundamental  spheroidal  form. 


Figure  3.7:  A  density  plot  of  the  diffraction  limited  intensity  psf,  h(v,u),  it 
is  rotationally  symmetric  around  the  u  axis. 

For  graphical  analysis  it  is  appropriate  to  look  at  three  different  views  or  cross- 
sections.  The  first  is  a  two-dimensional  spatial  psf  in  the  plane  perpendicular  to  the  optical 
axis  at  focus,  h(v,0).  Shown  in  Figure  3.8,  this  is  a  cross-section  of  an  infocus 
monochromatic  point  source.  It  displays  the  expected  amount  of  image  degradation  or 
blurring  applied  to  spatially  adjacent  coordinates  and  is  a  limiting  factor  in  determining  the 
spatial  resolution.  Substituting  u=0  into  equation  (3.8): 

'2/,(vr2 


h(v,  0)  = 


(3.13) 
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Figure  3.8:  A  theoretical  plot  of  h(v,0)  in  the  focal  plane. 

When  rotated  around  the  optical  axis  this  pattern  is  known  as  the  airy  disk.  The 
majority  of  the  intensity  falls  between  the  first  zeros  of  the  Bessel  function  at 

Ai-±1.22XF/#.  For  the  F/4  DOE  designed  for  the  prototype  Ar~5.7pm.  This  is  smaller 

than  the  1  Opm  pixels  of  the  CCD  detector.  Therefore,  once  detected,  the  infocus  diffraction 
limited  incoherent  psf  will  look  like  a  delta  function. 

The  second  view  of  interest  is  an  intensity  plot  through  focus,  h(0,u),  along  the 
optical  axis,  normally  considered  the  depth  of  focus.  Since  DOIS  interprets  z  as 
wavelength,  h(0,u.)  represents  the  expected  amount  of  blurring  applied  to  adjacent  spectral 
slices.  It  is  the  limiting  factor  that  determines  the  spectral  resolution  and  bandwidth. 
Substituting  v=0  into  equation  (3.9): 


h(0,  u)  = 


sin% 


(3.14) 


Figure  3.9:  A  theoretical  plot  of  h(0,u)  along  the  optical  axis. 
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In  Figure  3.9,  the  first  zero  in  intensity  along  the  axis  is  at  Az=±8XF/#2.  For  the 
prototype  Az  =  ±75fJm,  predicting  a  theoretical  incoherent  diffraction  limited  spectral 
resolution  of  AX^0.256/im. 


The  third  view  of  the  psf  is  along  the  geometrical  boundary,  where  u-v,  Calculated 
in  equation  (3.15)  and  plotted  in  Figure  3.10,  h(v,v)  dictates  the  blur  seen  from  coordinates 
that  are  different  in  both  spectra  and  spatial  location.  Substituting  u=v  into  equation  (3.9): 


Figure  3.10:  A  theoretical  plot  of  h(v,v)  where  u=v. 


(3.15) 


3.4.2  Experimentally  measured  incoherent  psf 

Unfortunately,  because  of  aberrations,  noise  and  other  errors,  the  blurring 
encountered  in  real  systems  is  greater  than  predicted  by  the  diffraction  limited  intensity 
point  spread  function.  Another  approach  to  obtain  the  three-dimensional  incoherent  psf  is 
by  direct  experimental  measurement  using  a  laser  illuminated  pinhole  to  simulate  a 
monochromatic  point  source.  The  following  discusses  the  validity  of  simulating  the 
incoherent  intensity  psf  with  images  of  a  coherent  pinhole. 

3.4.2. 1  Coherent  vs.  Incoherent  Imaging 

Coherent  imaging  is  described  by  the  amplitude  convolution  equation  [Goodman 

1968]: 
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*.•  -  ^  -  y„  Pobj  (xB,  y„  )dx0dyfl 


(3.16) 


where  Uim  represents  the  amplitude  distribution  in  image  space  from  the  amplitude 
distribution  in  object  space  Uobj,  and  A(x,,y,)  is  the  coherent  or  amplitude  impulse 

response.  However  it  is  actually  the  image's  intensity  Iim  which  is  detected: 

2 

jj %  -  x0,  y,  -  y„pohJ(x0, y0)dx0dy0  (3.17) 

Incoherent  imaging  is  a  linear  mapping  of  intensity  described  by  the  intensity 
convolution  integral: 


coherent:  4(w,)  = 


incoherent:  Iim(x„ y,)  =  JJ|%  -  xo, y,  - y0)\lobj{x0,  y0)dxBdyn 


(3.18) 


where  Iobj  is  the  object's  intensity  distribution  equal  to  |  Uoh]  | 2,  and  ^(v,,y;)2  is  the 

incoherent  or  intensity  impulse  response,  previously  represented  by  Mx^y). 

To  find  the  impulse  response,  the  input  object  U0bj(x0,y0)  is  a  perfect  point  source 

represented  by  a  delta  function,  S(xa,yJ.  Substituting  the  delta  function  into  equation 
(3.17)  for  coherent  imaging  yields: 

2 

4  (x„  y, )  =  Jj  %  -  x„,  y,  -  y„  )S(x0,  y„  )dx0dyo 
coherent:  ^3  j 

Likewise  for  incoherent  imaging,  substituting  the  delta  function  into  equation  (3.18) 
where  Iobj(x0,y0)  =  \5(x0,y0f  =  S(x0,y0): 
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incoherent: 


(3.20) 


2 

y,)  =  JJ|%  -  Vi  ~  y0)  %-  y0)dx0dy0 


=  %■>>',) 

Showing  that  if  the  input  object  was  indeed  a  delta  function  point  source,  the 
recorded  intensity  would  represent  the  incoherent  or  intensity  impulse  response, 


2 

h(xt,yt)=  h{xt,y)  ,  of  the  system.  However  in  practice  every  target  has  some  finite  width, 


thus  the  perfect  delta  function  point  source  does  not  exist. 

The  implications  of  the  difference  between  coherent  and  incoherent  illumination  can 
be  seen  if  the  object  being  imaged  is  a  pair  of  delta  function  point  sources  separated  by  the 
Rayleigh  distance  denoted  by  xr  [Gaskill  1978].  If  the  two  sources  are  incoherent,  the  total 
image  irradiance  is  simply  the  sum  of  the  image  irradiance  distributions  of  the  individual 
sources: 

incoherent:  Iim  (*, ,  y, )  = 

On  the  other  hand  if  the  two  sources  are  mutually  coherent,  the  overall  image 
irradiance  is  given  by  the  square  modulus  of  the  sum  of  the  amplitude  responses  of  the 


(3.21) 


individual  sources.  If  the  phase  difference  between  the  two  sources  is  <I>  than  the  irradiance 


may  be  expressed  as: 


%,y,) 

+  ejibh{xi  -  xT,  yt 

coherent:  = 
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+ 

h{xi  -xr,yt) 

(3.22) 


+  yi)h*{xi  -  xr,  y,)  +  ei<s>h*{xit  ?,-)£(*,■  -  y,) 


Notice  that  this  is  the  incoherent  result  of  equation  (3.21)  with  additional  interference 
terms.  It  is  these  interference  terms  which  create  the  difference  in  images  formed  from 
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coherent  and  incoherent  targets  and  the  difference  between  the  images  of  a  laser  illuminated 
pinhole  and  the  desired  incoherent  intensity  point  spread  function. 

Reynolds  et  al. ,  depict  the  difference  between  the  coherent  and  incoherent 
illumination  of  a  slit  object  in  the  book.  The  New  Physical  Optics  Notebook:  Tutorials  in 
Fourier  Optics,  which  is  reprinted  in  Figure  3.11.  The  theoretical  intensity  distribution  of 
the  coherent  and  incoherent  images  of  a  slit  is  shown  in  (a),  notice  the  coherent  illumination 
fringes  caused  by  the  interference  terms  of  equation  (3.22). 


Figure  3.11:  (a)  The  theoretical  intensity  distribution  in  the  coherent  and 
incoherent  images  of  a  slit,  (b)  the  recorded  image  of  a  coherent  slit, 
(reprinted  from  Reynolds  1989,  page  115) 

For  this  dissertation,  the  image  distribution  from  a  coherently  illuminated  5|im 
diameter  pinhole  is  detected,  which  is  a  circularly  symmetric  version  of  the  slit.  The 
diffraction  limited  intensity  distribution  in  the  coherent  and  incoherent  images  of  various 
pinholes  was  calculated  using  the  diffraction  limited  coherent  and  incoherent  impulse 
responses  from  Figure  3.12. 
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Figure  3.12:  The  diffraction  limited  impulse  response  for  (a)  coherent  and 
(b)  incoherent  illumination. 

Figure  3.13  shows  cross-sections  of  various  magnified  pinholes  in  the  left  column 
(MT=-.l),  alongside  the  resulting  diffraction  limited  intensity  distributions  from  coherent 
and  incoherent  illumination  plotted  together  in  the  right  column.  Notice  that  for  a  large 
160pm  pinhole,  Figure  3.13(a),  the  distribution  resembles  the  slit  image  of  Figure  3. 1 1(a). 
The  incoherent  image  (thin  line)  is  a  smoothed,  slightly  expanded  version  of  the  pinhole's 
rectangle  function,  where  as  the  coherent  image  (bold  line)  has  additional  ringing  about  the 
incoherent  distribution  due  to  the  interference  terms  of  equation  (3.22).  For  smaller 
pinholes  such  as  the  20pm  in  (d),  the  image  distribution  looks  less  like  the  slit  image  and 
more  like  the  incoherent  intensity  impulse  response.  For  the  5pm  pinhole  used  in  the 
experimental  psf  measurement  shown  in  (e),  the  intensity  distribution  matches  the  intensity 
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impulse  response  and  the  difference  between  the  coherent  and  incoherent  illumination  is 
negligible. 

Another  experimental  consideration  is  that  the  detector  pixels  are  1  Opm  wide.  This 
creates  a  scenario  where  all  of  the  distributions  of  Figure  3.13  are  unresolved.  The 
diffraction  limited  images  of  pinholes  smaller  than  100pm  fall  entirely  within  the 
boundaries  of  one  pixel.  Therefore,  the  interference  effects  of  the  coherent  source  aren't 
resolved  and  the  coherent  illumination  can  accurately  simulate  incoherent  illumination. 


3.4.2.2  Experimentally  measured  results 

A  set  of  images  of  the  coherent  illuminated  5pm  pinhole  taken  at  different  levels  of 
defocus  is  shown  in  Figure  3.14.  For  this  data  set  the  CCD  was  stepped  along  the  optical 
axis  with  a  step  size  of  dz=0.125  mm.  This  series  shows  the  effect  of  aberrations  on  the 
incoherent  point  spread  function.  The  predominate  aberration  is  spherical,  visualized  as 
lack  of  symmetry  about  best  focus.  The  spherical  aberration  caustic  creates  a  hard  ring  of 
intensity  on  the  marginal  side  of  focus  and  a  softer,  more  evenly  distributed  blurring  on  the 


paraxial  side  of  focus.  Theoretically  the  coherent  source  will  causes  interference  fringes  on 
top  of  the  aberrated  intensity  distributions,  however  again  they  are  unresolved  by  the  10pm 
pixels  and  average  to  an  appropriate  representation  of  the  incoherent  intensity  impulse 

response. 


marginal  <—  best  focus  — >  paraxial 


Figure  3.14:  Series  of  images  (32  x  32  pixels)  of  a  5  pm  pinhole 
illuminated  by  a  HeNe  laser  recorded  with  DOIS  through  focus. 

A  measured  3D  incoherent  intensity  point  spread  function  was  assembled  from 


stacking  32x32  pixel  images  of  a  pinhole  illuminated  with  a  542nm  GreNe  laser  at  32 
stages  of  defocus.  The  incremental  step  size  was  dz=0. 125mm,  corresponding  to  a 
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wavelength  step  size  dX=0.25nm.  A  cross-section  of  the  psf  cube,  h(r,z),  is  plotted  in 

Figure  3.15  as  (a)  a  density  plot  and  (b)  a  three-dimensional  plot.  The  3D  distribution  is 
circularly  symmetric  about  the  r=0  optical  axis. 


r  [mm] 


Figure  3.15:  The  experimentally  measured  intensity  point  spread  function, 
h(r,z),  (a)  a  density  plot  (b)  a  3D  plot. 


The  measured  spectral  intensity  point  spread  function,  h(z(X))=h(0,u)  ,  is  plotted  in 


Figure  3.16  with  the  theoretical  h(z(X))  determined  from  both  geometric  and  diffraction 

limited  theory.  The  z  axis  labels  include  the  spectral  conversion.  The  full  width  at  half  max 
(FWHM)  represents  the  spectral  bandwidth  of  DOIS.  The  FWHM  diffraction  limited 


bandwidth  is  AA=0.06nm  and  the  measured  spectral  bandwidth  is  AX=1.5nm,  both 


determined  at  X— 542nm.  See  Chapter  5  for  more  on  the  experimental  characterization  of 
DOIS. 
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Figure  3.16:  The  theoretical  and  measured  on-axis/spectral  intensity  psf, 

h(z(l)=h(0,u). 

The  three-dimensional  intensity  psf  of  Figure  3.15  can  be  Fourier  transformed  to 
yield  the  three-dimensional  optical  transfer  function  (OTF)  shown  in  Figure  3.17.  The  OTF 
is  required  for  several  of  the  digital  image  restoration  techniques  in  Chapter  8.  It  is 
important  to  note  that  there  is  a  significant  difference  in  the  behavior  of  the  experimental 
and  theoretical  OTFs.  Thus  it  is  extremely  important  to  use  the  experimental  OTF  wherever 
possible.  As  there  is  probably  less  difference  between  similar  DOE  lenses  than  there  is 
between  experimental  and  theoretical  functions,  it  is  reasonable  to  simply  interpolate  from 
experimental  curves  to  correct  for  changes  in  the  design  parameters. 
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Figure  3.17:  The  experimentally  measured  OTF(p,Q,  (a)  a  density  plot  (b)  a 
3D  plot. 
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4.  THREE  DIMENSIONAL  IMAGE  FORMATION  WITH  A 
DIFFRACTIVE  OPTICAL  ELEMENT  (DOE) 


The  Diffractive  Optic  Image  Spectrometer  (DOIS)  approach  to  image  spectrometry 
has  been  broken  into  the  series  of  coordinate  transformations  and  operators  that  are 
depicted  in  the  flow  chart  of  Figure  4.1,  [Sitter  1990,  Frieden  1967,  Barrett  1981]  The  top 
box  represents  imaging  with  the  diffractive  optical  element  (DOE)  in  DOIS  and  is  described 
in  section  4.1.  The  DOE  performs  two  operations,  Gaussian  imaging  with  a  Spectral 
Gaussian  Coordinate  Transform  and  image  degradation  characterized  by  Three-dimensional 
Shift  Variant  Transfer  Theory.  The  image  space  is  then  sampled  by  the  CCD  camera. 

The  acronyms  used  below  include  coordinate  transform  (CT),  shift  variant  (SV) 


and  shift  invariant  (SIV). 


object  =>  obj(x0 ,  y0 ,  l0 )  at  z0 
U-  Spectral  Gaussian  CT 
Gaussian  image  =>  o{x,  y,  z ) 

11  3D  SV  Transfer  Theory 
image  =>  i(xt ,  y, ,  z, )  _ 


<— narrow  waveband 


full  spectrum-* 


approximate"isotome  '  shift  invariant  imaging  shift  invariant  model  for  f  shift  variant  imaging 

,  ,  11  Isotome  CT  J 

11  SIV  Reconstruction  isotomic  image  =>  (  (*,•  ,  y,-  ,  z,- ) 

Gaussian  image  =>  o(x,  y,  z)  H  SIV  Reconstruction 

11  Inverse  Spectral  Gaussian  CT  isotomic  Gaussian  image  =>  o'  ( V ,  y' ,  z' ) 

reconstructed  object  =>  objr(x0,y0,  A0)  atz„  H  Inverse  Spectral  CT 

reconstructed  object  =>  objr(x0,y0,X0 )  atz0| 

Figure  4.1:  A  flow  chart  depicting  the  manipulation  of  the  three- 
dimensional  object  cube  as  it  is  imaged,  and  reconstructed. 

Image  processing  can  be  applied  to  reconstruct  the  object  under  one  of  two 

conditions,  shown  as  two  possible  paths.  Selection  of  a  path  will  depend  on  the  extent  of 


the  waveband  desired.  The  left  path  can  be  taken  within  a  narrowband  image  volume  called 
an  isotome  and  is  discussed  in  section  4.2.  In  an  isotome  the  point  spread  function  (psf)  is 
approximately  unchanging  and  is  considered  shift  invariant  for  modeling  and  reconstruction 
algorithms.  It  will  be  shown  in  section  4.2.2  that  achievable  isotomes  have  a  bandwidth  up 
to  20nm  in  the  visible  range  prototype  built  for  this  dissertation. 

When  the  spectra  of  interest  includes  the  entire  visible  range  it  is  necessary  to 
perform  the  coordinate  transform  (CT)  shown  in  the  right  path  of  Figure  4.1  and  discussed 
in  section  4.3.  The  CT  converts  the  image  space  with  a  shift  variant  (SV)  psf  into  an 
alternate  coordinate  system  with  a  shift  invariant  (SIV)  psf.  The  coordinate  transform 
designed  for  this  system  is  called  the  Isotome  Coordinate  Transform.  Once  applied, 
reconstruction  algorithms  using  linear  SIV  filtering  techniques  can  be  applied  to  obtain  an 
isotomic  Gaussian  image.  The  final  step  is  to  again  transform  the  data.  This  time  an  Inverse 
Spectral  Coordinate  Transform  is  used  to  reconstruct  the  object. 

4.1  Imaging  with  a  Diffractive  Optic 

Although  one  normally  thinks  of  a  diffractive  optical  element  (DOE)  as  a 
monochromatic  imaging  device,  they  can  be  used  to  image  polychromatic  sources.  Sweatt 
has  shown  a  mathematical  equivalence  between  a  diffractive  optical  element  and  an  ultra 
high  index  lens  [Sweatt  1979],  The  DOE  performs  as  a  lens  with  severe  chromatic 
aberration.  It  is  this  aberration  which  provides  the  dispersion  for  DOIS. 

It  is  important  to  remember  that  all  imaging  systems  are  imperfect,  introducing 
aberrations  which  distort  the  image.  Although  high-quality,  diffractive/refractive  optical 
kinoforms  can  be  designed  to  minimize  image  distortion  and  aberrations  such  as  spherical 
and  coma  [Stone  1988],  the  DOIS  spectral  sectioning  technique  is  plagued  with  defocus 
from  the  chromatic  aberration.  This  results  in  errors  in  the  recorded  spectral  slices  which 
have  to  be  corrected  with  post-detection  image  processing. 


60 


Figure  4.2:  A  diffractive  optical  element  imaging  various  spectral  point 
sources. 

A  good  place  to  start  this  discussion  is  to  consider  imaging  a  point  source.  Figure 
4.2  shows  a  DOE  imaging  point  sources  of  various  wavelengths.  Notice  that  the  DOE  does 
several  things.  The  spectra  is  spread  out  along  the  optical  axis  and  the  magnification  is 
different  in  each  focal  plane.  These  functions  are  described  by  Gaussian  imaging. 
Additionally,  the  images  of  the  point  sources  are  not  perfect  delta  functions.  Each  is  a 
spheroidal  intensity  distribution  around  the  "ideal"  Gaussian  image  point.  This  blurring  is 
described  by  three-dimensional  transfer  theory.  The  spheroid  intensity  distributions 
represent  the  three-dimensional  point  spread  functions  of  the  DOE,  and  cause  both  spatial 
(transverse)  and  spectral  (longitudinal)  blurring. 

4.1.1  Spectral  Gaussian  Coordinate  Transform 

The  imaging  can  be  represented  by  a  coordinate  transform,  performing  both  a 
scaling  function  and  transforming  the  spectral  data  into  a  spatial  coordinate,  z.  Starting  with 
the  imaging  equation: 

4--UI  (4.1) 

So  Z«  Z 
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where  fB  is  the  focal  length  at  Xa,  z„  is  the  object  distance  and  z  is  the  Gaussian  image 
distance.  Introducing  the  diffractive  optic  spectral  response  [Swanson  1989]: 

fo^o  ~  fill'd  ~  D 

fo=£y!L  =  T  (4.2) 

x„  K 

D  is  the  DOE  design  constant  where  fd  is  the  design  focal  length  at  the  design  wavelength, 
Xd.  Substituting  (4.2)  into  (4.1)  results  in  the  diffractive  optic  imaging  equations: 


i  =  ±  +  I 

D  zn  z 


X  _  A  £±*« 


<=>  z  =  - -  z„ 

l  X„zn-D) 


The  conversion  between  spectra  and  image  distance,  z,  at  a  fixed  object  distance,  z0,  is 
described  by  equation  (4.4). 

The  transverse  and  longitudinal  magnification,  MT  and  ML,  are: 


M-r  = 


z„  Xnzn  -  D 


and  M,  =  -  —  =  —M. 


Notice  that  the  magnification  is  dependent  on  the  design  constant,  D,  and  the  fixed  object 
distance,  z0,  leaving  wavelength,  as  the  only  variable. 

The  DOE  transforms  the  object  objr(x0,y0,X0)  at  z0,  into  a  Gaussian  image, 
o(x,y,z),  with  the  following  Spectral  Gaussian  Coordinate  Transform’. 

o{x,y,z)  =  obj[MTx0,MTy0,-MTz0)  (4.6) 


f  ~D  1  f  —D  D 

o(x,y,z)  =  obj  - - -  - - -  y0,  — - -  z0 

Xnzn -D  I  \  X?  -D  Xnzn  - D 
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4.1.2  Three-dimensional  Transfer  Theory 

As  previously  mentioned,  the  actual  image  of  a  point  source  formed  by  a  DOE  is  a 
blurred  version  of  the  Gaussian  point,  seen  in  Figure  4.2.  The  blurring  is  represented  by 
Three-dimensional  Transfer  Theory,  based  upon  the  assumptions  of  (1)  superposition  and 
(2)  stationarity.  It  is  well  known  that  these  conditions  are  met  in  two-dimensions  by  an 
object  which  radiates  incoherently  from  an  "isoplanatic"  area  of  an  object  plane.  These 
conditions  can  be  extended  to  three-dimensional  imaging  within  "isotomic"  volumes  of  the 
object  space  [Frieden  1967].  Frieden's  theory  is  extended  below  to  describe  the  DOIS 
three-dimensional  spectral/spatial  system. 

The  theory  of  superposition  implies  that  since  the  object  is  assumed  to  be 
incoherent,  the  intensities  from  different  elements  of  the  object,  o(x,y,z),  are  additive.  The 

total  intensity  in  the  image,  i'(x,,y(,Z,),  is  given  by: 

oo  oo  oo 

i{xjt  y„  Z.HJJ  h(x,  y,  z;xn  y,.,  z,)o(x,  y,  z)dxdydz  (4.8) 

where  h{x,y,z\x ^y^z)  is  the  space  variant  incoherent  impulse  response  or  psf.  This 

impulse  response  physically  represents  the  intensity  distribution  in  image  space  from  the 
illumination  by  a  monochromatic  point  source  in  object  space.  It  is  particularly  important 
for  estimating  the  image  degradation  caused  by  surrounding  coordinates.  The  three- 
dimensional  psf  predicts  cross-talk  and  blurring  from  coordinates  both  within  the  same 
spectral  image  and  from  surrounding  spectra. 

Classical  two-dimensional  transfer  theory  defines  an  "isoplanatic"  plane  where  the 
psf  is  slowly  varying  and  approximately  stationary.  For  DOIS  an  "isotomic"  volume  of  the 
object  or  Gaussian  image  space,  called  an  isotome,  can  be  defined  where  the  three- 
dimensional  psf  is  dependent  on  the  displacement  from  the  Gaussian  point  and  not  the 
location  of  the  Gaussian  point  itself,  described  by  equation  (4.9). 
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(4.9) 


h(x,y,z;Xi,yi,Zi)  =  h(xi-x,yi-y,zi-z) 

Physically,  this  requires  that  the  psf  must  remain  invariant  under  changes  in 
position  of  the  Gaussian  point.  This  will  be  true  if  each  point  within  an  isotome  has 
essentially  the  same  wave  aberrations  and  magnification. 

Image  space  can  be  broken  into  several  isotomic  volumes,  V.  Each  isotome,  which 
represents  a  narrow  waveband,  can  be  treated  by  linear  shift  invariant  filtering  theory.  If  a 
particular  spectral  image  or  a  narrowband  spectral  curve  is  needed,  the  detection  and 
processing  can  be  confined  to  an  isotome  following  the  left  path  of  Figure  4.1. 

4.2  Isotomic  Imaging  in  a  narrow  waveband 

Substituting  equation  (4.9)  into  (4.8),  the  way  in  which  an  arbitrary  incoherent 
spectral  image  is  distorted  can  be  predicted  from  the  resulting  intensity  convolution 
equation.  It  describes  the  convolution  of  the  Gaussian  image  with  the  incoherent  psf  in  an 
isotome. 

oo  oo  oo 

Kxi » >/»*/)=  J  J  J  hixi  -  x’  yt  ~  y>  z,  ~  z)o(x,  y,  z)dxdydz  (4.10) 

in  short  hand  notation: 

Kxi ’  .  z« )  =  h(x,  y,  z)  *  *  *  o(x,  y,z )  (4.11) 

Applying  the  Fourier  Integrals: 

0(£  V>Q  =  in 

I{t¥.  f)  =  j  (4.12) 

H&yr.C)  =  J  JJ  h(x,y,z)e“^”t(‘)dMlydz 
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The  Fourier  transform  of  the  psf,  H(  £  y/,  Q,  is  called  the  optical  transfer  function 
(OTF),  or  in  three-dimensional  microscopy  texts,  the  contrast  transfer  function  (CTF). 

l^,xi/,C)  =  H^,y/,C)0(^C)  (4.13) 

Here  capitals  are  used  to  refer  to  the  Fourier  transforms  and  the  convolution  is  now 
simplified  to  a  multiplication.  In  matrix  form,  dropping  the  variables: 

I  =  HO  (4.14) 

If  the  isotome  object  and  psf  are  known,  the  image  can  be  modeled  with  equation  (4.14), 
finishing  with: 

i  =  InverseFourier{I} 
i(xi ,  y,  ,Z,)  =  InverseF ourier  {/(£,  lf,Q} 

This  will  simulate  the  image  without  the  errors  and  noise  introduced  by  detection. 

4.2.1  Reconstruction 

Within  the  isotome,  if  the  recorded  spectral  images,  i(xi,yi,zi),  the  object  distance, 
z0  and  the  OTF  are  known,  an  Inverse  Shift  Invariant  Filter  and  Inverse  Spectral  Gaussian 
Coordinate  Transform  can  be  applied  to  reconstruct  the  object,  objr(xo,y0,f0)  at  za.  The 
above  transfer  theory  can  be  applied  in  reverse  using: 


0  =  H-1I 

(4.16) 

o  =  InverseFourier{0} 

(4.17) 

(;c,;y,z)  =  InverseFourier^O[t,,\j/,lf)\ 

where  I=/( £  y/,  Q  is  the  Fourier  transform  of  the  image,  i(xi,yj,zi),  and  H'1  is  the  inverse 

of  the  OTF  matrix  which  when  applied  as  an  inverse  filter  performs  the  required 
deconvolution. 
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To  convert  o(x,y,z)  into  a  reconstructed  spectral/spatial  object,  objr(x0,y0,X0)  at 
z0,  the  following  Inverse  Spectral  Gaussian  Coordinate  Transform  has  been  derived. 

(4.18) 

obJr{x0,y0,K)  =  o  -~x,-—y,D  (4.19) 

<  Z  Z  V  ZZo  )  j 


0bjr(Xo^„X)  =  0\~,^-,D\-  +  - 
{  Mt  Mt  U 


4.2.2  Define  isotome  bandwidth 


To  determine  which  path  of  Figure  4.1  should  be  taken,  it  must  be  determined  if  the 
image  cube,  or  bandwidth  of  the  image  cube,  can  be  modeled  as  a  shift  invariant  isotome. 
To  define  an  isotome  the  invariance  on  h  must  be  quantified  [Frieden  1967],  Refer  to 
section  3.4.1  for  the  theoretical  description  of  the  incoherent  point  spread  function  h..  The 
invariance  over  the  isotome  can  be  measured  by  the  invariance  of  h(u,o)  and  h(0,v)  under 

changes  Az  or  AX.  This  can  be  quantified  by  the  differential  radial  distance  Ar=Rv  to  the 

first  zero  of  h(v,0 )  in  equation  (3.13),  and  the  axial  distance  Az=  Zu  to  the  first  zero  of 
h(0,u)  in  equation  (3.14): 

Ar=R,  =±\.22l{ffj  or  (4.20) 

*=Z‘=±K£)  °X=±2-f-  (4'21) 

recall  that  D  is  the  design  constant  fdXd.  The  differentials  ARv  and  AZU  measure  the  variance: 


ARV  =  Az,X 

AZU  _  A z  -t  _  A M 

Zu  zt  M 


(4.22) 

(4.23) 
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In  practice,  Frieden  defines  an  isotome  as  a  region  where  the  magnification  changes 
less  than  10%  (AML  /ML<10%).  Therefore,  over  small  spectral  regions  or  small  amounts  of 

axial  translation  the  image  space  can  be  considered  isotomatic  and  SIV  linear  filtering 
algorithms  will  hold.  Using  the  differentials  in  equation  (4.22)  and  (4.23),  which  are 
plotted  in  Figure  4.3,  see  that  for  a  target  at  an  image  distance  of  zo>10f  and  a  bandwidth, 

AAA,  up  to  5%,  that  the  variance  AR/R<1%  and  AZ/Z<6%.  This  indicates  that  the 
conditions  for  an  isotome  hold  within  a  spectral  region  up  to  a  20nm  bandwidth. 


4.3  Full  Spectrum  Imaging 

The  technique  of  section  4.2  can  be  applied  to  individual  isotomes  to  model  and 
reconstruct  data  from  a  narrowband  source.  However,  when  the  spectra  of  interest  includes 
the  entire  visible  range,  it  is  necessary  to  perform  the  coordinate  transform  shown  in  the 
right  path  of  Figure  4.1,  which  creates  a  shift  invariant  model  for  shift  variant  imaging  . 

Barrett  and  Sitter  introduced  a  technique  to  model  the  shift  variant  blurring  on  three- 
dimensional  images  [Barrett  1981,  Sitter  1990].  The  model  consists  of  a  coordinate 
transformation  of  the  object  followed  by  a  shift  invariant  blur  operation,  which  is  in  turn 
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followed  by  a  second  coordinate  transformation  of  the  image.  This  technique  can  be 
rederived  to  apply  to  the  reconstruction  of  the  object  when  the  image  space  is  shift  variant. 

Recall  from  equation  (4.8),  if  an  imaging  system  is  not  isotometric,  it  is  modeled  by 
the  shift  variant  integral  of  the  form: 

*(*/>  ?,<*,)=  JJ  J  h(x>  y> z; xi’ yp  zi )°(x>  y>  z)dxdydz 


where  i(xi,yi,zi)  represents  the  image  space  distribution,  o(x,y,z)  is  the  Gaussian  object 

distribution  and  is  the  space  variant  incoherent  psf.  Sitter  showed  that  for 

incoherent  imaging  systems,  there  exist  coordinate  transformations  from  object  space  and 
image  space  to  an  intermediate  space  (which  will  be  denoted  by  primed  quantities)  where 
the  integral  of  equation  (4.8)  can  be  expressed  as  a  convolution.  That  is,  there  exists  a 

coordinate  transformation  from  i(*„y,  ,z,  )  to  i1  (x/ ,y/ ,z/)  and  from  o(x,y,z)  to 
o'  (x'  ,y',z ')  so  that: 


*'(*/'. ’>z,')  =  J  J \h'(x',y',z';xi',yi',zi')o'(x’,y,,z')dx'dy,dz'  (4.24) 


where:  H  (x'  ,y'  ,z'  ;*/,?/  ,z/)  =  H  (*/ ,y/  -y'  ,z/  -z) 

resulting  in: 


(4.25) 


*'(*,- yt Zi ')  =  J  J  J h’(x,  '-x',  y.  '-y',  Zi  '-z')o'(x',  y',  z')dx’dy'dz'  (4.26) 


where  i  [x-  ,y-  ,z-)  and  o'  (x'  ,y'  ,z')  are  the  coordinate  transformed  irradiance  distributions 
and  h'  is  the  incoherent  psf  in  the  transformed  space. 


4.3.1  Isotome  Coordinate  Transform 

A  coordinate  transform  was  derived  for  the  DOE  to  transform  the  shift  variant 
image  space  i  to  the  shift  invariant  isotome  V  with  the  assumption  that  aberrations  are 
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minimized  and  that  magnification  is  the  major  cause  of  the  space  variance.  The  isotome 
coordinate  transform  applies  an  inverse  magnification  transforming  a  polygon  shaped  space 
with  variable  magnification  into  a  unit  magnification  space  pictured  in  Figure  4.4. 


Figure  4.4:  Visual  representation  of  the  isotome  coordinate  transform 
forming  an  image  space  with  unit  magnification 

In  order  to  remove  the  variable  magnification  which  is  a  function  of  z/A),  the 


following  coordinate  transform  is  used: 


Mt  =  — — 


Ml  =  - 


( ^ 

\ZoJ 


=  —Mi 


*i  y, 

^  Mr  ’  Mt'  Mh 


t{x;,y:,z;)=i\[ 


(4.27) 


(4.28) 


h(x',y',z')  =  i 


(4.29) 


Equations  4.28  and  4.29  completely  reverse  the  magnification  of  the  Gaussian 
imaging.  However  in  some  applications  it  may  be  practical  to  transform  the  image  cube  to  a 
intermediate  magnification,  such  as  to  match  the  magnification  of  the  measured  image  cube 
to  the  magnification  of  a  measured  point  spread  function. 
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4.3.2  Reconstruction 


Now  that  the  image  space  is  essentially  all  one  isotome,  the  entire  full  spectrum 
object  cube  can  be  reconstructed  with  the  inverse  transfer  reconstruction  that  was  used  in 


section  4.2.1,  substituting  the  primed  variables: 

i'{xi y, Zi ')  =  h'(x',  y',  z')***  o'(x',  y',  z ')  (4.30) 

r  {£,ij/,C)  =  If  (Z,wX)0'  (4.31) 

O^H-’r  (4.32) 

o'  =  InverseFourier{  O’ } 

o'  (x'  ,y',z')  =  InverseF ourier  (£,(//,£)}  ^ 

4.3.3  Inverse  Spectral  Coordinate  Transform 


The  magnification  from  the  xt  and  yt  coordinates  was  removed  in  the  coordinate 
transform  of  equation  (4.28).  The  final  coordinate  transform  reconstructs  the  spectral 
information  Xn  from  z-,  using  equation  (4.4). 


at  z 


o 


°bjr(Xo’  y<n  K)  =  o' 


(4.34) 
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5.  PROTOTYPE  CHARACTERIZATION 


translation 

stage 


lib  jack 


5.1  DOIS  Prototype 

A  prototype  of  the  DOIS  described  in  Chapter  2  was  assembled  and  characterized  at 
Rome  Laboratory’s  Photonics  Center.  Figure  5.1  is  a  photo  of  the  DOIS  prototype.  The 
DOE  is  on  the  right  and  connected  to  the  CCD  camera  with  a  tubular  baffle.  Although  a 
system  designed  for  field  use  would  include  electronic  control  translation  stages,  the 
laboratory  prototype  is  manually  controlled.  The  DOE  is  mounted  on  an  Onel  rail  earner 
which  slides  on  a  60  cm  Oriel  rail  marked  with  a  metric  scale  for  location  measurement. 


Figure  5.1:  Photo  of  the  DOIS  prototype. 

The  camera  is  first  secured  to  a  micropositioning  translation  stage  for  fine  control, 

then  on  another  rail  carrier.  Each  carrier  is  marked  with  a  graduated  rule  for  positioning  on 
the  rail  within  1  mm,  and  the  micropositioning  vernier  on  the  detector  is  accurate  to  .001 
inches  (25.7  |Xm).  Manual  translation  of  the  camera  or  DOE  is  possible.  With  a  target  at 
infinity,  the  lens  would  be  moved  without  adversely  impacting  the  object  distance  s0  and 
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the  detector  to  DOE  distance  sa  for  each  wavelength  would  be  the  spectral  focal  distance, 

six  =  However,  for  target  distances  restrained  by  the  dimensions  of  the  laboratory,  the 

best  results  are  obtained  by  keeping  the  DOE  stationary  and  moving  the  detector.  The 
detector  carrier  is  stepped  with  coarse  increments,  dz  =  1  mm,  corresponding  to  spectral 

step  size  dA  ranging  from  lnm  to  3nm,  at  400nm  and  650nm  respectively.  For  finer 

spectral  measurements,  the  detector  carrier  was  placed  at  a  set  position  and  the 

micropositioning  vernier  was  tuned  with  dz  =  .005”,  corresponding  to  a  dA  of  .lnm  to 

.4nm  depending  on  the  wavelength  region. 

5.2  The  Diffractive  Optic  Element 

The  DOE  fabricated  for  this  prototype  is  a  single  etch  Fresnel  Zone  Phase  Plate. 
Figure  5.2  lists  the  figures  of  merit  for  the  DOE  in  DOIS.  The  dominant  effect  of  only  a 
single  etch  DOE  is  the  low  first  order  diffraction  efficiency.  This  limits  the  amount  of  light 
which  is  properly  imaged  and  increases  the  minimum  detectable  target  radiance,  but  it 
performs  sufficiently  for  this  proof  of  concept,  laboratory  demonstration.  A  well  designed 

DOIS  would  incorporate  a  high  efficiency  phase  grating  with  r\  ~  100%  at  the  design 

wavelength.  The  DOE  could  be  a  multilevel,  N  >  16  photolithography  phase  plate,  or  a 
continuous  profile  grating  generated  by  holography,  diamond  turning  or  a  laserwriter  as 
discussed  in  Chapter  3. 
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DOE  type 

2  level,  phase  etched,  zone  plate 

Fabrication  at  the  National 
Nanofabrication  Facility  of  Cornell 
University 

ebeam  mask  generation 
contact  lithography  exposure 
ion  reactive  etched 

substrate 

BK-7  optical  flat,  3mm  thick 

diameter 

2”,  5  cm  clear  aperture 

design  focal  length  f0 

200  mm 

design  wavelength  A,0 

588  nm 

theoretical  diffraction  efficiency  ,  r\ 

40  %  at  588  nm 

experimentally  measured  rj 

27%  at  670  nm  laser  diode 

35  %  at  632.8  nm  HeNe 

30  %  at  542  nm  GreNe 

MTF  @  Ve 

32  lp/mm 
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The  ring  pattern  of  Figure  5.3  was  recorded  by  illuminating  DOIS  with  a  high 
intensity  HeNe  beam.  The  first  order  focus  is  super  saturated  and  is  seen  on  the  bottom 
left.  The  DOE’s  zero  order  diffraction  efficiency  is  approximately  5%  at  633nm.  The 
resulting  zero  order  beam  has  no  optical  power,  passing  straight  through  the  lens.  It  adds 
background  photon  noise  but  is  normally  not  intense  enough  to  be  detected  as  signal. 
However,  when  illuminated  with  the  laser  beam,  the  undiffracted  light  projects  a  shadow  of 
the  zone  plate  on  the  detector. 


The  first  order  diffraction  efficiency  is  a  measure  of  the  percentage  of  the  light 
collected  that  is  imaged  to  the  desired  detector  location.  This  figure  of  merit  dictates  the 
minimum  target  radiance  that  will  be  detected.  The  setup  in  Figure  5.4  was  used  to 
experimentally  determine  the  DOE  s  diffraction  efficiency.  To  determine  the  spectral 
dependence  of  the  diffraction  efficiency  the  DOE  was  characterized  at  three  wavelengths.  A 
laser  diode  was  used  which  emitted  at  670  nm  and  measurements  were  obtained  at  632.8 
nm  and  542nm  with  two  Helium  Neon  lasers.  Figure  5.5  shows  the  experimental  results 
plotted  along  with  the  theoretical  spectral  first  order  diffraction  efficiency  [Chapter  3]. 

The  laser  beam  was  directed  into  an  integrating  sphere/detector  apparatus  and  the 
power  was  recorded,  PtotaI  Iaser.  The  DOE  was  then  placed  in  the  beam  and  the  integrating 
sphere  was  aligned  to  collect  only  the  first  diffracted  order  and  again  the  power  was 
recorded,  Pistorder.  The  diffraction  efficiency  is  the  ratio  of  these  two  values. 


p 

—  1 st  order 

'Ist  order  n 

1  total  laser 


(5.1) 


To  account  for  the  efficiency  at  different  radii  and  zone  spacing,  r\  was  determined  with  the 


beam  hitting  the  DOE  at  half  radius  and  repeated  at  almost  full  aperture.  The  two  radii 
measurements  were  averaged. 
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Figure  5.4:  Diffraction  efficiency  experimental  setup  (a)  determine  the  total 
power,  (b)  first  order  power  at  r/2.  (c)  first  order  power  at  full  aperture. 
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The  DOE  MTF  was  measured  at  the  Diffractive  Optics  Test  Facility  at  Rochester 
Photonics  Corporation,  RPC.  This  is  an  automated  testbed  that  measures  the  MTF  of  a 
diffractive  optic  using  a  double  knife  edge  technique  [Faklis  1993].  Figure  5.6  (a)  shows 
the  measured  tangential  and  sagital  MTFs  of  the  DOE  at  540  nm,  as  well  as  the  diffraction 
limited  MTF.  The  dramatic  drop  from  the  DC  term  is  due  to  the  low  diffraction  efficiency 
of  the  DOE.  The  second  curve  (b)  is  the  MTF  normalized  with  the  diffraction  efficiency. 
Warren  Smith,  in  his  book  on  lens  design  [Smith  1992],  gives  the  goal  of  a  MTF  of  50%  at 
50  lp/mm  for  high  quality  photographic  lenses  and  20%  at  30  lp/mm  over  90%  of  the  field. 
This  DOE  falls  well  below  this  requirement,  however  Figure  5.6(c)  is  the  measured  MTF 
of  a  well  designed,  high  efficiency  DOE  lens  and  is  presented  to  show  that  this  technology 
can  greatly  exceed  standard  design  requirements.  Additionally,  the  DOE  design  can  be 
optimized  using  computer  design  codes  such  as  ZEMAX  to  reduce  aberrations  both  on  and 
off  axis.  The  substrate  surface  curvatures  are  used  for  minimizing  aberrations  by  producing 
a  kinoform,  which  adds  a  low  power  refractive  surface  for  balancing  aberrations  such  as 
spherical  and  coma. 
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1.1- 


Figure  5.6:  DOE  MTF, 
measured  at  RPC  Diffractive 
Optic  Testing  Facility. 
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(a)  measured  and  diffraction  limited 
MTF  at  540  nm. 


(b)  Normalized  546  nm  MTF, 
MTF/r|. 


(c)  MTF  performance  of  a  well 
designed  high  efficiency  DOE 
measured  at  559  nm,  provided  as 
an  example  of  the  technology’s 
capability. 


5.3  The  CCD  camera 


The  CCD  camera  in  the  prototype  is  a  SONY  XC-75  whose  specifications  are 
found  in  Figure  5.7. 


SONY  CCD  Monochrome  Video  Camera 

XC-75 

Pickup  device 

Interline-transfer  CCD 

effective  picture  elements  (pixels) 

768  x  494  (horizontal  /  vertical) 

sensing  area 

1/2  inch 

CCD  drive  frequency 

15.7  kHz  vertical 

14.3  kHz  horizontal 

pixel  size 

8.4  x  9.8  pm  (horizontal  /  vertical) 

chip  size 

7.95  x  6.45  mm  (horizontal  /  vertical) 

scanning 

525  lines,  2: 1  interlace 

gain  settings 

AGC  automatic  gain  control 

F  fixed  or  M  manual 

video  S/N  ratio 

56  dB 

minimum  illumination 

3.0  lux 

sensitivity 

400  lux  with  AGC  ON 

Figure  5.7:  DOIS  camera  specifications. 


The  CCD  camera  was  used  without  the  conventional  lens  attachment,  and  the  DOE 
was  mounted  «  200  mm  away  from  the  camera  with  a  homemade  baffle  of  concentric 
tubes,  depicted  in  Figure  5.1,  blocking  most  stray  photons  from  the  room.  The  output  of 
the  camera  was  connected  to  a  monitor,  a  framegrabber  board  and  a  video  recorder.  The 
images  were  recorded  digitally,  one  image  at  a  time.  The  computer  has  the  capability  of 
recording  many  frames  consecutively,  but  memory  constraints  prevented  this.  The  low 
efficiency  of  the  DOE  made  it  necessary  to  have  the  automatic  gain  control  on  for  most  of 
the  experimental  demonstrations.  This  is  unfortunate  because  the  capability  of  absolute 
intensity  readings  is  lost.  However,  the  AGC  did  a  wonderful  job  of  automatically 

adjusting  to  the  spectral  efficiency  T|^  of  the  DOE  and  camera.  In  an  improved  system,  the 
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gain  will  be  set  constant  and  the  intensity  will  be  corrected  by  dividing  the  spectral  intensity 
data  by  the  DOE  spectral  efficiency ,  t\x,  and  the  camera  spectral  response  curve,  measured 

below. 

The  camera  spectral  response  was  measured  using  a  Tungsten  Halogen  lamp  and  an 
Oriel  grating  monochrometer.  Figure  5.8  is  the  relative  spectral  response  of  the  CCD 
camera  in  its  original  packaging.  This  curve  suggested  that  there  was  a  filter  on  the  camera, 
possibly  to  mimic  the  response  of  the  human  eye.  After  careful  removal  of  the  faceplate,  a 
filter  was  found  between  the  protective  window  and  the  CCD  chip.  Figure  5.9  is  the 
spectral  response  of  the  same  CCD  with  the  unknown  filter  removed.  The  image  of  Chapter 
2  were  recorded  with  this  filter,  but  the  analytical  demonstration  in  the  next  chapter  was 
performed  after  the  filter  was  removed. 


spectral  curve  of  the  lamp  and  monochrometer  grating  divided  out. 
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Camera  Response  without  filter 


Figure  5.9:  the  experimentally  determined  spectral  response  curve  of  the 
SONY  CCD  after  the  removal  of  a  factory  placed  filter.  The  dashed  line  is 
the  raw  data,  and  the  solid  line  has  the  spectral  curve  of  the  lamp  and 
monochrometer  grating  divided  out. 


5.4  DOIS  prototype  characteristics 

Since  imaging  is  a  large  part  of  the  requirements  of  the  system,  and  many  readers 
may  not  have  seen  the  imaging  capabilities  of  a  simple  Frensel  Zone  Plate,  a  256  x  256 
image  of  the  Air  Force  Resolution  Target  is  included  as  Figure  5.10.  The  resolution  target 
used  was  made  of  chrome  on  glass.  It  was  illuminated  with  the  Tungsten  Halogen  lamp 
and  filtered  with  a  10  nm  line  width  filter  at  550  nm.  The  resolution  target  was  placed  over 
ground  glass  so  that  the  target  and  not  the  filament  was  imaged.  The  resolution  target  was 
12.5  x  12.5  mm,  S0=219.5  cm  and  S;  =  23.5  cm.  These  distances  correspond  to  a 
wavelength  of  552  nm.  The  x  and  y  extents  of  the  image  were  measured  in  MacPhase  to  be 
134  x  134  pixels,  corresponding  to  a  transverse  magnification  of  0.107  with  10|im  pixels. 


M„ 


=  _A  =  _i  =  _A 

so  x0  y„ 

-  _  23-  5cm  _  134 pixels  *  1 0| im 
219.5cm  12.5  mm 


-0.107 


(5.2) 
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These  measurements  confirm  that  the  spectral  imaging  and  transverse  magnification 
equations  hold  true  for  this  diffractive  optic. 


Figure  5.10:  A  pictorial  representation  of  the  system's  MTF,  a  DOIS  image 
of  an  Air  Force  Resolution  Target. 


DOE  offset 
-1.5  cm 


Detector  offset 
+  3  cm 


CCD 


iiiiiiiiiQiiiiiiiiiimiiiiiimiiimiiiiiiiiiimiimiiiimiiiiimiigiiiiimmmmiiiiii 

0 - - ruled  railing  i  ;  60  cm 


i|  DOE  scale 

L.  _ 


'^i~ ZDet'ZDOE 


Det  scale 
i 

- ► 


—  Figure  5.11:  Schematic  of  DOIS. 

It  is  necessary  to  determine  the  exact  location  of  both  the  DOE  and  detector  to 
properly  determine  the  spectral  value.  Figure  5.11  is  a  schematic  diagram  of  DOIS  and 

illustrates  the  distances  used  in  the  excel  spreadsheet,  Figure  5.12,  which  calculates  X  from 


S0,  the  DOE  and  detector  positions.  Precise  offsets  between  the  carrier  markings  and  the 
DOE  and  CCD  surfaces  are  required.  These  were  first  roughly  determined  from  the 
mounting  parameters  and  then  calibrated  with  the  known  HeNe  and  GreNe  wavelengths. 
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Once  set  at  -1.5  cm  and  +3  cm,  they  remained  accurate  for  a  variety  of  source  targets, 
wavelengths  and  distances. 
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Figure  5.12:  Excel  spreadsheet  calibrating  the  DOE  and  detector  carrier 
locations  z,  to  a  corresponding  wavelength,  X. 


5.5  Spatial  and  Spectral  point  spread  function  (psf) 

Figure  5.13  is  a  series  of  images  (32  x  32  pixels)  of  a  5  pm  pinhole  illuminated  by 
a  HeNe  laser  recorded  with  DOIS  through  focus.  The  separation  between  each  image  is 

dz=  .005”,  which  is  approximately  a  spectral  step  size  of  dA,=  0.4  nm.  Each  image  is 

plotted  with  its  own  min,  max  plot  range  to  maximize  the  contrast  and  boldly  show  the  spot 
diagram  pattern.  The  geometrical  image  at  best  focus  is  contained  within  one  pixel,  so  this 
will  experimentally  represent  the  system's  spatial  psf  for  various  stages  of  defocus.  These 
images  not  only  tell  us  the  amount  of  blur  added  to  our  image,  but  they  also  graphically 
demonstrate  the  dominant  lens  aberration,  third  order  spherical.  The  spot  diagrams 
correspond  to  under  corrected  spherical  aberration.  Notice  the  soft  even  blur  as  the  images 
progress  into  focus,  but  then  the  hard  ring  of  light  after  focus.  Figure  5.14  is  a  collection  of 
the  cross  sections  of  each  of  these  psfs  grouped  as  “approaching  focus”  and  “leaving 
focus”  respectively,  and  plotted  as  greyscale  intensity  value  verses  pixel  number.  The 
spectral  images  are  each  32x32  pixels.  Notice  that  the  255  greylevel  maximum  of  the  psf  is 
centered  on  pixel  number  16  and  falls  of  quickly  to  a  greylevel  of  zero  within  6  pixels  on 
either  side  (shown  are  pixel  numbers  10  through  22,  beyond  them  the  psf  is  zero). 
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Figure  5.13:  Series  of  images  (32  x  32  pixels)  of  a  5  pm  pinhole 
illuminated  by  a  HeNe  laser  recorded  with  DOIS  through  focus. 
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Figure  5.14:  Spatial  Impulse  response  h( r)  at  632.8  nm,  plotted  as 
greylevel  vs.  pixel  number  for  locations  (a)  approaching  focus,  (b)  leaving 
focus. 

Another  important  system  parameter  can  be  obtained  from  this  data  set.  Plotting  the 
intensity  distribution  of  one  pixel  through  focus  becomes  a  measure  of  the  influence  of  one 
slice  of  the  data  on  the  next  plane  of  data.  In  3D  microscopy  this  measures  the  degradation 
in  image  quality  of  one  spatial  plane  z  on  the  next,  but  in  this  configuration  it  is  a  measure 
of  the  degradation  of  one  wavelength  by  the  surrounding  spectral  channels.  This  translates 
into  the  spectral  impulse  function  which  determines  spectral  bandwidth  or  resolution  of  the 

system.  The  experimentally  determined  coherent  spectral  impulse  function  h(k)  is  plotted  in 

Figure  5.15,  the  intensity  distribution  of  the  center  pixel  through  focus  at  (a)  632.8  nm 
HeNe  and  (b)  542  nm  GreNe. 
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Figure  5.15:  Spectral  Impulse  response  h(k),  plotted  as  intensity  of  one 
pixel  through  focus  (a)  for  632.8  nm.  (b)  for  542  nm. 
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Figure  5.16:  The  experimentally  determined  spectral  resolution. 

These  spectral  psfs  can  be  scaled  with  wavelength  to  translate  through  the  entire 

visible  range.  Each  step  size  is  dz=.005”=.127  mm,  HeNe  dz/zx =.0 1 27/20.93=. 06%, 
GreNe  dz/zx=. 0127/24.91=. 05%.  To  use  the  experimentally  determined  GreNe  psf  data 

increments  of  dz/zx=.05%  should  be  used.  For  example,  to  filter  properly  at  450  nm, 
z=29.4  cm,  the  step  size  dz  should  be  equal  to  0.15  mm  ~  .006”. 
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6.  IMAGE  FORMATION 


The  previous  chapters  have  concentrated  on  describing  DOIS's  operation 
theoretically.  At  this  point  it  is  useful  to  try  and  verify  that  the  theory  does  indeed  predict 
the  image  formation  function  of  DOIS.  This  Chapter  presents  both  computer  generated 
(CG)  and  experimentally  determined  image  sets  of  known  targets.  As  you  will  see,  the 
linear  system  theory  of  Chapter  4  accurately  predicts  output  images  formed  by  DOIS.  The 
CG  images  are  very  similar  to  those  experimentally  measured,  which  verifies  that  the 
theory  is  appropriate. 

6.1  Mercury  577/579  nm  doublet  with  an  X  aperture 

An  excellent  example  of  DOIS’s  performance  is  its  ability  to  resolve  a  two 
nanometer  doublet  (577/579nm)  from  a  Mercury  Lamp.  To  test  both  spectral  and  spatial 
imaging,  a  known  aperture  was  placed  in  front  of  the  lamp.  The  aperture  was  two  crossed 
slits  in  metal  that  make  an  X.  The  lamp  was  positioned  at  the  same  object  distance  that  the 
GreNe  pinhole  psf  images  were  recorded  at  and  images  were  recorded  every  dz=0. 125mm 

corresponding  to  a  spectral  step  size  of  dX=0.3nm. 

The  algorithm  in  section  4.2  was  implemented  in  Mathematica  and  used  to  simulate 
the  experimental  scenario.  Figure  6. 1  highlights  the  Mathematica  code  for  this  model.  It 
creates  a  three-dimensional  (32x32x32)  image  matrix,  from  2  three-dimensional  input 
matrices,  a  spectral  object  matrix  and  the  measured  psf  matrix  [Figure  2.12].  Each  is 


formed  by  stacking  32  images  that  are  32  by  32  pixels  in  dimension. 

OB  J=centeredFourier[Flatten[object]]  (6.1) 

OTF=centeredFourier[Flatten[psf]]  (6.2) 

IMAGE  =  OTF  *  OBJ  (6.3) 

image=Partition  [  Abs  [centeredln  ver  seFourier  [IMAGE]  ]  ] _ (6.4) 


Figure  6.1:  Mathematica  code  that  implements  the  linear  system  algorithm. 
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For  this  simulation  it  was  assumed  that  Mercury  emits  two  narrow  lines  at  577  nm 
and  579  nm.  Pictured  in  Figure  6.2  is  the  series  of  assumed  input  objects,  paired  with  the 
CG  output  images  that  were  generated  from  the  Mathematica  linear  system  algorithm,  along 
side  the  experimentally  measured  images.  The  input  object  was  assumed  to  have  non-zero 
values  in  channels  7,8,9  and  14,15,16  exclusively,  corresponding  to  wavelengths  of 
576.6,  576.9,  577.2,  578.7,  579.0  and  579.3  nm.  The  corresponding  object  spectral 
radiance  curve  is  Figure  6.3(a)  showing  emission  only  around  577  nm  and  579  nm.  The 
experimentally  measured  OTF  matrix  from  Figure  2.14  was  used. 

The  images  created  by  this  algorithm  accurately  depict  the  image  formation 
properties  of  DOIS.  In  Figure  6.2  you  can  see  the  excellent  correlation  between  the  CG  and 
the  experimentally  recorded  images  of  the  third  columns.  The  spectral  radiance  curve  of  the 
CG  images  is  Figure  6.3(b).  As  expected  the  psf  does  degrade  the  spectral  separation  of  the 
doublet.  Yet,  the  doublet  can  still  be  detected  and  can  be  restored  with  various 
deconvolution  techniques  [Chapter  8],  The  corresponding  spectral  curve  from  the 
experimental  image  set  is  Figure  6.3(c),  which  shows  even  better  spectral  separation. 
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Figure  6.2:  Simulated  spectral  object  &  image  channels  of  a  Mercury 
loublet  X,  using  experimentally  measured  psf(x,y,z(A-)). 


Figure  6.3:  Spectral  radiance  plots  of  a  Mercury  doublet  577/579nm  (a)  the 
assumed  input  source  object,  (b)  the  computer  generated  results,  (c)  the 
experimental  result. 


6.2  Spectral/spatial  analysis 


The  above  target  demonstrated  high  spectral  resolution  by  being  able  to  discern  a 
2nm  doublet  while  imaging  the  narrow  slits  of  the  X.  Small  targets  encounter  only  minimal 
blurring;  however  as  the  target  increases,  so  does  the  spectral/spatial  blurring.  This  is 
investigated  by  generating  CG  images  of  targets  of  varying  spatial  extent  with  the  Fourier 
model  using  the  measured  OTF. 


The  data  can  be  displayed  in  many  ways.  A  conventional  representation  is  a  series 
of  spectral  images  like  Figure  6.2,  or  a  spectral  curve  at  one  pixel,  Figure  6.3.  However  to 
visualize  the  spectral/spatial  relationship,  the  following  discussion  will  utilize  a 
spectral/spatial  cross-section  image.  Figure  6.4  shows  two  three-dimensional  data  cubes. 


The  cross-section  is  formed  as  a  zxx  or  zxy  plane  , shown  as  the  images  above  the  cubes. 


Notice  that  the  monochromatic  zxy  image  looks  like  a  slit  and  the  image  from  the 

doublet  emission  looks  like  a  double  slit.  The  height  of  the  slit  represents  a  dimension  of 
the  aperture  of  the  target  and  the  width  is  the  spectral  bandwidth  of  the  emission  line. 
Figure  6.5  presents  the  spectral/spatial  on-axis  cross-sections  of  the  point  spread  functions 
used  during  this  project.  The  contrast  has  been  enhanced  to  show  features. 
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Figure  6.4:  3D  object  cube  and  zy  vie’1 
(a)  monochromaticly  (b)  a  doublet  soi 


Figure  6.5:  Cross-sections  of  the  three-dimensional  point  spread  functions 
used  in  CG  image  formation  and  object  reconstruction. 
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6.2.1  Monochromatic  target,  increasing  aperture 

Linear  system  theory  was  used  to  model  the  images  formed  by  DOIS  of  a 
monochromatic  target  with  six  rectangular  apertures  of  different  sizes.  The  cross-sections 
are  shown  in  Figure  6.6.  The  images  were  generated  by  forming  a  3D  matrix,  32x32x32 
where  the  only  non-zero  elements  are  a  square  of  ones  in  the  16th  plane,  simulating  a  square 
monochromatic  target.  The  size  of  the  square  is  increased  from  1  pixel  to  the  full  32x32 
plane  (1,  2, 4,  8,  16  and  32  pixels  respectively).  As  expected,  as  the  target  size  is  increased 
its  total  radiance  is  increased,  causing  increased  blurring  in  both  the  spectral  and  spatial 

directions.  Each  two-dimensional  image  is  a  32x32  zxy  cross-section  of  the  total  three- 

dimensional  cube  with  the  corresponding  spectral  curve  plotted  below  it. 


Figure  6.6:  Computer  generated  monochromatic  target  of  increasing  size. 
6.2.2  Doublet  target,  increasing  aperture 


A  more  interesting  target  is  the  same  increasing  square  aperture  in  front  of  a 
doublet.  The  images  in  Figure  6.7(a)  are  the  CG  assumed  input  sources  with  the  square 
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aperture  at  z=14  and  20.  The  series  in  Figure  6.7(b)  are  the  CG  images  which  were 
calculated  with  the  measured  OTF.  For  comparison  a  mercury  lamp  with  577nm  and  579 
nm  line  emissions  illuminated  various  apertures.  The  measured  z,y  images  were  extracted 
from  the  measured  cubes  and  shown  in  Figure  6.7(c).  Notice  the  strong  agreement  between 


the  simulated  and  experimental  images  and  curves. 
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Looking  at  the  experimentally  measured  image  of  the  50pm  and  100pm  apertures 
on  the  left  side  of  Figure  6.7,  notice  that  the  target  drifted  to  a  lower  pixel  in  the  y  direction 
as  the  detector  was  scanned  along  the  z  axis.  This  misregistration  visualizes  the  effect  of 
both  a  changing  magnification  and  off  axis  aberrations  which  skew  the  orientation  of  the 
point  spread  function,  Figure  6.8.  This  artifact  can  be  eliminated  with  a  constant 
magnification  design  which  is  optimized  for  minimum  aberrations. 


? _ 10 _ 20 _ 90 _ 40 _ 50 _ 60 

spectral  step  #  2 


Figure  6.8:  A  point  spread  function  askew  from  off-axis  aberrations. 

Another  instance  to  be  discussed  while  investigating  targets  of  increasing  size  is 
when  the  target  is  image  to  the  edge  pixels  of  the  CCD.  If  the  target  is  imaged  to  the  edge  of 
the  detector  array  field-of-view  (FOV),  the  information  about  the  target  which  is  blurred 
outside  of  the  FOV  will  be  lost.  Additionally,  the  simulation  and  reconstruction  algorithms 
are  based  on  transfer  theory  assuming  a  shift-invariant  point  spread  function. 
Mathematically  this  assumes  that  the  psf  and  OTF  matrices  are  circulant,  that  each  row  has 
equivalent  elements  where  the  values  are  shifted  and  the  end  value  wrap  around  to  the  other 
side  as  shown  in  equation  (6.5). 


94 


(6.5) 


abed 
d  a  b  c 
c  d  a  b 
b  c  d  a 


Figure  6.9  shows  the  wrap  around  in  the  y  direction  that  is  assumed  if  the  psf  is  on  the 
edge  of  the  field.  A  field  stop  is  proposed  to  limit  the  FOV  to  be  smaller  than  the  useable 
detector  area.  Pixels  outside  of  the  FOV  will  record  the  blurred  intensity  from  targets 
imaged  to  the  edge  of  the  FOV.  By  recording  the  fully  blurred  target  the  reconstruction  will 


be  improved. 


D _ 10 _ 2D _ 90 _ 40 _ 50 _ 6D  ^ 

spectral  step  #  z 

Figure  6.9:  The  psf  of  a  pixel  at  the  edge  of  the  detection  area. 
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7.  PROTOTYPE  DEMONSTRATION  WITH  FOUR  TARGETS 


The  DOIS  Prototype  was  demonstrated  with  four  known  multi-spectral/spatial 
targets.  The  sources  include  a  Tungsten  Halogen  lamp,  a  Mercury  lamp,  a  HeNe  laser  and 
a  GreNe  laser.  For  complete  control,  each  source  illuminated  a  known  template  which 
became  the  target.  Figure  7.1  is  a  view  of  the  targets  imaged  by  a  conventional  lens  and 
recorded  with  the  monochrome  CCD  camera. 


Tungsten  Halogen 

Mercury 

1  x  1  mm  eross 
with  10  gm  slits 

HeNe 

50  |am  pinhole 

v/ 

5  |am  pinhole 

■m 

- 

'V 

GreNe 

10  jam  pinhole 

Figure  7.1:  The  four  targets:  a  Tungsten  Halogen  Lamp,  a  Mercury  Lamp,  a 
HeNe  and  a  GreNe. 

7.1  Testbed  Specifications 

The  Tungsten  Halogen  lamp  is  a  broadband  source,  as  shown  in  Figure  7.2.  It 
illuminated  a  50  |im  pinhole,  creating  a  several  pixel  object  and  preventing  images  of  the 
filament.  The  spectral  emission  of  the  Mercury  lamp  is  shown  in  Figure  7.3.  Mercury 
makes  an  excellent  spectrometer  target;  it  is  a  multi-line  emitter  with  six  major  lines  in  the 
DOIS  spectral  range.  Of  particular  interest  are  a  pair  of  emission  lines,  called  a  doublet,  at 
577nm  and  579nm.  DOIS  can  resolve  this  doublet,  demonstrating  that  the  system’s  spectral 


resolution  is  better  than  2nm!  The  aperture  used  with  the  Mercury  lamp  is  a  1  x  1  mm 
cross,  x,  with  10  pm  slits.  The  remaining  sources  were  two  Helium  Neon  lasers.  The 
HeNe  emits  at  632.8  nm  and  the  GreNe  at  542  nm.  They  illuminated  a  5^im  and  10pm 
pinhole  respectively,  simulating  monochromatic  point  sources. 


relative 

spectral 


radiance 


reylevel 


400  500 

wavelength  [nm] 


(a) 


_ wavelength  [nm] 


(b)  step  size:  dz  =  0.1cm,  dA,  =  0.8nm-03nm 

Figure  7.3:  Spectral  radiance  curve  of  the  Mercury  Lamp  (a)  reprinted  from 
Applied  Optics,  (b)  measured  with  DOIS. 
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The  four  targets  were  placed  an  equal  distance  from  DOIS  and  beamsplitters  were 
aligned  to  create  the  illusion  that  they  were  all  in  the  same  object  plane.  Figure  7.3  is  the 
schematic  of  the  demonstration  testbed.  The  Tungsten  Halogen  lamp  was  directly  in  line 
with  DOIS.  The  Mercury  and  laser  sources  were  reflected  into  the  field  of  view  with 
beamsplitters.  The  object  distance  s0  was  91.5”  and  the  DOE  was  held  constant.  The  CCD 
camera  was  manually  scanned  and  the  detector  carrier  location  and  inch  scale  readings  were 
recorded  for  each  spectral  image.  The  256x256  monochromatic  images  of  each  target  at 
different  wavelengths  is  given  in  Figure  7.5. 


Figure  7.4:  Schematic  of  the  four  target  testbed. 
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365nm  Mercu 

Q 

'  383nm  Tungsten  Halogen  Lamr 

)  404nm  Lamp  &  Mercu 

4 1 5nm  Tungsten  Halogen  Lamp  435nm  Lamp  &  Mercury  470nm  Tungsten  Halogen  Lam 


542nm  Lamp,  Mercury  &  GreNe  546nm  Lamp  &  Mercury  560nm  Tungsten  Halogen  Lam 


578  Lamp  &  Mercury  600nm  Tungsten  Halogen  Lamp  632.8nm  Lamp  &  HeNe 

Figure  7.5:  Sample  of  the  key  spectral  DOIS  images  in  this  demonstration. 
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7.2  Results 

DOIS  was  used  to  image  the  multi  spectral/spatial  targets.  As  previously  discussed, 
the  DOE  was  kept  stationary  and  the  detector  was  stepped  through  a  wide  range  of 
distances.  As  the  detector  was  moved,  the  various  targets  popped  in  and  out  of  the  image 
depending  on  their  spectral  content,  as  seen  in  Figure  7.5.  However,  since  Tungsten 
Halogen  is  a  greybody  emitting  over  the  entire  visible  range,  the  image  of  the  50pm  pinhole 

was  present  continuously  in  the  z  range  for  A,  from  383nm  to  650nm.  The  Tungsten 

pinhole's  intensity  grew  at  longer  wavelengths  and  its  center  shifted  several  pixels  because 
of  the  magnification  changes.  The  best  focus  of  each  spectral  line  was  determined  and  the 
corresponding  detector  locations  were  recorded.  The  table  in  Figure  7.6  shows  the 
calculation  of  the  emission  line  peak  wavelengths  from  this  experiment.  The  results 
correspond  perfectly  to  the  known  emissions  of  the  lasers  and  the  Mercury  lamp  of  Figure 
7.4,  and  are  remarkably  precise  considering  the  manual  stepping  and  reading  of  the  carrier 


location. 


source 

s0  [cm] 

DOE 

carrier 

Det  carrier 
[cm] 

Inch  scale 
[in] 

Si  [cm] 

emission  line 
peak  wavelength 
[nm] 

HeNe 

232.41 

7 

24 

0.508 

20.2 

633 

Hg  Lamp 

232.41 

7 

26.1 

0.515 

22.3 

579 

Hg  Lamp 

232.41 

7 

26.1 

0.481 

22.4 

577 

Hg  Lamp 

232.41 

7 

27.5 

0.4875 

23.8 

546 

GreNe 

232.41 

7 

27.7 

0.504 

23.9 

542 

Hg  Lamp 

232.41 

7 

34.3 

0.488 

30.6 

435 

Hg  Lamp 

232.41 

7 

36.8 

0.51 

33.0 

407 

Hg  Lamp 

232.41 

7 

36.8 

0.396 

33.3 

404 

Hg  Lamp 

232.41 

7 

41 

0.43 

37.4 

365 

Figure  7.6:  Peak  spectral  lines  of  the  target,  calculated  from  the  locations 
of  best  focus. 
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Attempts  were  made  to  experimentally  verify  the  Mercury  lamp’s  spectral  emission 
curve  with  the  grating  monochrometer.  Of  particular  interest  is  the  doublet  at  577  nm  and 
579  nm.  The  signal  from  the  lamp  was  low  and  the  smallest  useable  exit  slit  width  was  «  .2 
mm  corresponding  to  a  bandwidth  of  5  nm  making  the  results  unusable.  However,  this 
emphasizes  the  usefulness  of  even  this  inefficiently  designed  DOIS. 

To  determine  the  spectral  radiance  curves  of  the  targets,  series  of  images  were 
recorded  with  known  increments  of  dz.  The  spectral  radiance  is  defined  as  the  intensity 
value  of  the  image  at  various  wavelengths.  A  C  program  was  written  to  open  each  image 
and  record  the  value  of  a  specific  pixel  in  a  list.  This  list  was  plotted  against  wavelength 
creating  a  measured  spectral  radiance  curve.  A  different  approach  needed  to  be  implemented 
when  the  change  in  magnification  throughout  the  spectral  range  caused  the  image  to  drift 
across  several  pixels.  Under  this  condition,  a  window  was  defined  around  the  image 
location  and  the  value  extracted  for  each  image  is  the  maximum  pixel  value  in  the  window. 

The  experiment  was  executed  twice,  each  time  with  a  different  spectral  resolution. 
The  first  data  run  was  low  resolution,  testing  the  entire  visible  spectrum.  It  incorporated 
manual  translation  of  the  detector  carrier  in  increments  of  dz=0.1mm.  The  inch  scale 
vernier  was  at  a  constant  setting.  The  second  was  high  resolution,  concentrating  on  a 
limited  spectral  range  around  each  emission  line.  The  detector  carrier  was  kept  stationary 
and  the  inch  vernier  was  adjusted  in  0.005”  or  0.01”  increments. 

The  results  from  the  dz=0.1mm,  low  resolution  recording  of  the  Mercury  Lamp  are 
plotted  in  Figure  7.4(b).  There  is  excellent  agreement  between  the  measured  and  known 
spectral  emission  lines  of  the  characteristic  curve  7.4(a)  from  Applied  Optics[ Levi  1968]. 
However,  the  DOE  diffraction  efficiency  [Figure  5.5]  weighted  the  relative  intensity  of 
each  spectral  line. 

The  results  from  the  Tungsten  Halogen  pinhole  are  found  in  Figure  7.7.  Included 
are  6  spectral  images  clipped  from  the  full  images  to  show  just  the  50pm  pinhole.  The 
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spectral  radiance  plot  is  constructed  from  data  recorded  in  low  resolution  steps.  The 
measured  data  points  are  plotted  as  *  and  a  line  represents  the  best  fit.  An  artifact  of  the 
changing  magnification  is  that  the  image  drifts  from  a  condition  of  being  centered  on  a  pixel 
to  one  where  the  image  straddles  several  pixels,  as  shown  in  Figure  7.8.  This  alignment 


problem  explains  the  ripple  effect  on  the  intensity  verses  wavelength  curve. 


X  =  400nm  450nm  500nm  550nm  600nm  650nm 

(a) 

200  F  Tungsten  Halogen  — ^ — . 


greylevel 


*  .*  * 


saturation 


data  points 
best  fit  to  data 


400  450  500  550  600  650 

wavelength  [nm] 

_ (b) _ 

Figure  7.7:  Tungsten  Halogen  Data  (a)  spectral  images  of  the  50pm 
pinhole,  (b)  the  measured  spectral  radiance  curve. 
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Another  difficulty  of  this  data  run  was  that  the  signal  beyond  525nm  saturated  the 
detector  array.  A  more  accurate  plot  of  the  Tungsten  Halogen  is  found  in  Figure  7.3(b). 
The  experiment  was  run  again  with  a  larger  source  whose  image  was  several  pixels, 
avoiding  pixel  hopping,  and  with  the  gain  set  to  prevent  saturation.  The  results  accurately 
measure  the  spectral  radiance  when  compared  to  the  handbook  curve  in  Figure  7.3(a). 

The  results  from  the  GreNe  and  HeNe  laser  sources  are  shown  in  Figures  7.9  and 
7.10.  In  addition  to  the  spectral  radiance  plots,  a  representative  series  of  images  is  given  in 
each  figure.  The  images  are  32x32  windows  which  have  been  clipped  from  the  256x256 
recorded  images.  As  expected  DOIS  recorded  narrowband  high  intensity  spikes  for  both. 
Figure  7.9  shows  the  HeNe  narrowband  spectral  curve  centered  at  632.9nm.  Figure  7.10 
shows  the  same  for  the  GreNe  at  542nm.  Notice  that  with  constant  spatial  steps, 
dz=0.005”,  the  change  in  magnification  causes  the  spectral  step  size  to  change  within  the 

spectral  range.  At  634nm,  the  spectral  steps  were  0.4nm,  yet  at  541.5nm,  dk=0.2nm. 


measured  spectral  radiance  curve  (dz=0.005”=0.0125cm,  d?i=0.3-0.4nm). 
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measured  spectral  radiance  curve  (dz=0.005”=0.0125cm,  d?i=0.2-0.3nm). 


The  Mercury  target  was  examined  a  second  time  with  small  step  sizes  around  each 
Hg  emission  line.  Figure  7.11  shows  the  resulting  spectral  curves  for  the  365,  404,  435, 
546  and  577/579nm  lines  with  the  clipped  infocus  spectral  images. 

Close  observation  of  the  576.9nm  image  in  Figure  7.1 1(e)  reveals  a  focused  cross 
with  a  highly  defocused  cross,  and  a  the  center  that  looks  like  a  diamond.  Moving  to 
579nm,  the  diamond  comes  to  focus  as  the  center  of  the  579nm  image,  and  the  577  cross 
defocuses  and  disappears.  This  corresponds  to  the  computer  generated  images  in  Chapter 
6.1.  This  recorded  data  around  the  577nm  and  579nm  Hg  doublet  shows  that  the  doublet  is 
resolved  without  processing.  However,  still  greater  discrimination  is  obtained  with  the 
processing  applied  in  Chapter  8. 


105 


350  360  370  380 

_ wavelength  [nm] 


size:  dz-0.01”=0.Q25cm.  diU.2nm 


b)  step  size:  dz=0.01”=0.025cm  dX=0.3nm 


_ wavelength  |nm] _ 

(d)  step  size:  dz=0.005”=0.0125cm.  dk=0.3nm 


T 


X  =576.9 


577.2  577.5 

greylevel 

200  r 


577.8  578.1 

Hg  577/579im 


579.0  nm 


560  570  580  590  600 

wavelength  [nm] 

„.  (e)  steP size:  dz=,005  d?v=0.3nm  "  ' 

figure  7.11:  The  measured  spectral  radiance  curves  of  each  Mercury  line  with  the  in-focus 
32x32  ptxe  spectral  image(s)  of  the  center  of  the  cross,  (a)  365nm  line,  (b)  404nm  line, 
(c)  435nm  line,  (d)  546nm  line,  (e)  577/579nm  doublet. 
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8.  OBJECT  RECONSTRUCTION  ALGORITHMS 


The  results  in  Chapter  7  show  that  DOIS  performs  image  spectrometry.  The  image 
quality  and  spectral  resolution  available  directly  from  DOIS  will  be  useful  in  many 
applications.  Notwithstanding,  applying  an  object  reconstruction  algorithm  can  improve 
DOIS’s  performance  by  removing  the  blur  from  surrounding  channels,  resulting  in  a  more 
accurate  representation  of  the  target's  object  cube.  Several  object  reconstruction/ 
deconvolution  techniques  are  presented  here.  They  have  been  applied  to  both  computer 
generated  (CG)  and  experimentally  measured  3D  DOIS  image  cubes. 

In  Chapters  1  and  4,  DOIS  was  described  as  a  diffractive  spectral  sectioning  device. 
The  output  of  the  DOE,  described  by  the  Spectral  Gaussian  Coordinate  Transform,  is  an 
image  space  which  is  mathematically  equivalent  to  the  output  image  space  of  a  three- 
dimensional  imaging  microscope,  described  as  spatial  sectioning  microscopy. 

There  are  proven  digital  image  restoration  techniques  used  in  3D  microscopy  that 
eliminate  blurring  from  out-of-focus  adjacent  slices  of  the  3D  object.  It  follows  that  these 
same  techniques  can  be  adapted  and  applied  to  DOIS  images  to  remove  blurring  caused  by 
the  out-of-focus  spectral  images.  This  chapter  reviews  the  algorithms  which  have  been 
adapted  and  applied  to  DOIS. 

Since  not  all  applications  require  the  same  resolution,  three  general  processing 
techniques  will  be  discussed  in  this  chapter;  nearest  neighbor,  inverse  filtering  and 
constrained  iterative  deconvolution.  Each  approach  provides  different  amounts  of 
deblurring  at  various  computational  expense 

The  nearest  neighbor  algorithm  reconstructs  an  image  by  looking  only  at  the  impact 
from  the  adjacent  images  on  either  side.  This  requires  very  little  prior  knowledge  and  could 
be  applied  in  real  time. 
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Inverse  filtering  tries  to  eliminate  the  effect  of  the  instrument  function  on  the  output 
image  cube.  The  system  OTF  must  be  known  from  either  theoretical  prediction  or  direct 
measurement.  The  main  task  is  to  invert  the  OTF  and  multiply  the  image  set  by  this  inverse. 
Variations  include  applying  an  apodization  or  regularized  inverse  filter  for  noise 
suppression.  Additionally,  tools  such  as  singular  value  decomposition  may  be  required  to 
invert  the  large,  possibly  singular,  OTF  matrix.  Inverse  filtering  reconstruction  is  limited 
by  the  cutoff  frequencies  of  the  OTF;  frequencies  beyond  the  cutoff  can't  be  recovered. 

Iterative  techniques  which  apply  physical  constraints,  such  as  positivity,  to  the 
reconstructed  object  can  theoretically  reconstruct  frequencies  beyond  the  OTF  cutoff.  An 
additional  advantage  of  iterative  algorithms  is  that  they  don't  require  the  possibly  difficult 
task  of  inverting  the  OTF  matrix. 

These  algorithms  should  provide  equal  or  greater  reconstruction  to  the  DOIS 
spectral  sectioning  system  than  available  with  spatial  sectioning  3D  microscopy  because  of 
two  advantages.  First,  when  imaging  a  3D  spatial  object  such  as  a  photo-luminescent  cell, 
the  cell  itself  will  cause  blurring  and  information  loss  due  to  absorption  and  scattering 
through  the  object.  The  spectral  dimension  in  DOIS  won't  cause  this  degradation. 
Secondly,  3D  microscopy  has  the  demanding  task  of  removing  the  blur  from  slices  with 
very  low  z-spatial  resolution.  The  spatial  image  in  the  z  direction  is  normally  on  the  order 
of  the  psf  blur,  requiring  reconstruction  in  the  very  low  and  noisy  regions  of  the  OTF. 
However  in  DOIS  the  dispersion  is  so  great  that  the  total  spectral  range  translates  into  a  dz 
spatial  region  which  is  many  orders  of  magnitude  greater  than  the  spatial  z  width  of  the  psf, 
with  a  ratio  similar  to  the  radial  x,y  spatial  imaging.  Therefore,  in  most  situations 
reconstruction  around  and  beyond  the  OTF  cutoff  isn't  necessary. 

It  is  important  to  note  that  even  the  simplest  nearest  neighbor  method  produces  a 
substantial  improvement.  However  the  most  accurate  results  are  obtained  with  the  SVD 
Inverse  Filtering  algorithm. 
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8.1  Digital  Representation 

In  Chapter  4  the  object  and  image  cubes  were  represented  as  continuous  functions 
of  x,y,z  and  X.  However  as  seen  in  Figure  2.1,  the  images  are  sampled  in  x  and  y  by  the 

detector  array  and  in  z  by  stepping  the  array  along  the  optical  axis.  This  is  best  modeled  by 
discrete  object  and  image  functions.  In  practice  the  object,  image  and  psf  cubes  are  more 
accurately  represented  by  three-dimensional  matrices.  Recall  equation  (4.10): 

oo  oo  oo 

Kxi  >  y, ,  z, )  =  J  J  J  °ix .  *  zHxi  ~ x’  yi  -  *  z< ~  z)dxdydz 

this  equation  can  be  written  as  summations  rather  than  integrals  with  discrete  axis  intervals 
Ax,  Ay  and  Az. 

i(lAxi,pAyi,jAzi')  =  ^^^°{mAx,qAy,kAz)h(lAxi  -mAx,pAy(. -qAy,jAzt  —  A:Az)(8.1) 

m= 1  <?=1  k= 1 

using  subscripts  to  describe  the  Az  location: 

ij  (x,  y)  =  X  X  X  (x>  y)  *  *  *hi-k  (x’  y)  (8  -2) 

m= 1  q- 1  fc=l 

where  jc=  lAxt  and  mAx,  and  y=  pAyt  and  qAy.  The  value  of  j-k  represents  the  number  of 
steps  of  defocus. 

8.2  Nearest  Neighbor  Reconstruction 

The  Nearest  Neighbor  technique  is  routinely  used  in  three-dimensional  microscopy 
for  deblurring  optically  sectioned  image  sets  using  data  from  a  single  focal  plane  above  and 
a  single  focal  plane  below  to  correct  the  central  plane.  This  is  a  reasonable  approximation 
for  two  reasons.  First,  since  spectral  images  from  adjacent  focal  planes  will  contribute  most 
strongly  to  the  central  plane,  their  effects  are  most  important  to  consider.  Secondly,  these 
adjacent  image  planes  are  actually  a  collection  of  all  of  the  object  channels  with  various 


109 


amounts  of  defocus.  This  rationale  was  first  used  by  Castleman  [1979]  and  later  by  Agard 
[1989]  to  justify  using  only  the  nearest  neighbors: 

Simplifying  equation  (8.2)  by  dropping  the  x  and  y  variables,  models  the  system  as 
a  series  of  spatial  slices: 


N  N-j 

I \°k*hhk=  %ok+j*hk  (8.3) 

*=>  k=i-j 

This  simply  states  that  the  j  image  is  a  sum  of  convolutions  of  the  various  object  spectra 
with  the  appropriate  defocus  psfs.  Pulling  out  and  solving  for  the  in-focus  k= 0  term: 


ij=°j*ho+'iok+j*hk+  ^ok+j*hk 

k=\-j  *=1 


-1  N-j 

°j  *ho=ij  -  2>*+,  "2X;  *K  (8.5) 

k=l-j  *=1 

The  h0  is  the  in-focus  psf  of  the  DOE.  This  equation  states  that  the  object  at  j,  convolved 
with  the  in-focus  psf,  is  given  by  the  image  at  level  j  minus  a  sum  of  adjacent  object 
channels  that  have  been  blurred  by  out-of-focus  psfs  hk. 

This  suggests  that  the  object  at  channel  j  can  be  recovered  by  subtracting  from  the 
image  at  channel;  a  series  of  adjacent  objects  blurred  by  the  defocus  transfer  function.  If  a 
simultaneous  solution  approach  is  abandoned,  the  adjacent  object  channels  ob  .  are  not 
available.  However,  the  adjacent  images  ik+j  are  available.  From  equation  (8.4)  each  image 
contains  the  corresponding  object  channel  plus  a  sum  of  defocused  adjacent  object 
channels.  Furthermore  ignore  the  effect  of  the  in-focus  psf,  and  recognize  that  both  the 
recorded  image  j+I  and  image  j_,  contain  the  out-of-focus  contributions  from  the  entire 
stack.  In  the  positive  direction,  blurring  image  j+I  by  one  step  of  defocus  provides  a  good 
approximation  to  the  out-of-focus  contributions  that  contaminate  ij.  Combining  this  with  a 
similar  defocus  in  the  negative  direction: 
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(B.6) 


Oj~ij-c{ij-x*h-X+ij+x*K) 
where  h_,  and  h+]  are  psfs  that  approximate  the  blurring  due  to  defocus  by  one  step  in  either 
direction  (-Az  and  +Az),  and  c  is  an  adjustable  constant  (c=0.45  was  found  to  be 
appropriate).  This  suggests  that  the  defocused  structures  can  be  partially  removed  by 
subtracting  images  from  adjacent  planes  that  have  been  convolved  with  the  appropriate 
defocus  psf. 

The  nearest  neighbor  algorithm  of  equation  (8.6)  was  implemented  in  Mathematica 
and  applied  it  to  twenty  2D  spectral  images  of  the  Mercury  doublet.  The  code  is  shown  in 
Figure  8.1.  The  resulting  reconstructed  objects  are  shown  along  side  the  recorded  images 
in  Figure  8.2.  This  technique  significantly  reduced  the  image  blurring  and  greatly  improved 
the  spectral  resolution.  The  original  and  reconstructed  spectral  radiance  curves  of  a  pixel  are 
in  Figure  8.3.  The  dip  between  577nm  and  579nm  was  increased  from  20%  in  (a)  to  55% 
in  (b).  Additionally,  since  the  psf  is  non-symmetric  around  focus,  the  recorded  data  had  an 
incorrect  higher  value  at  579  nm.  The  use  of  separate  psfs  for  opposite  defocus  steps  Az— 
+1  and  -1  reconstructed  the  equal  intensities  of  each  line.  The  best  results  were  obtained 
with  the  section  3.4.2  experimentally  measured  psf  images,  one  from  each  side  of  the 
infocus  plane,  for  h  ,  and  h+1. 

The  nearest  neighbor  algorithm  is  particularly  useful  when  one  specific  spectral 
image  is  needed.  It  is  very  fast  and  requires  little  memory.  Computations  could  be  done 
using  generally  available  image  processing  hardware.  Thus,  it  should  be  possible  to  do  the 
all  of  the  required  calculations  in  essentially  real  time.  While  one  cannot  expect  this 
technique  to  recover  the  object  function  exactly,  it  does  promise  to  improve  images  and 
spectral  resolution  at  reasonable  expense. 
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For[z=l,  z<=20,  z++, 

IMAGE[[z]]=cfft[image[[z]]];] 

For[i=l,  i<=32,  i++, 

For[j=l,  j<=32,  j++, 

IMAGEPl[[I,ij]]=IMAGE[[l,ij]]OTFl[[ij]] 

IMAGEPl[[2,i,j]]=IMAGE[[2,i,j]]OTFl[[i,j]] 

IMAGEPl[[3,ij]]=IMAGE[[3,i,j]]OTFl[[i,j]]: 


IMAGEPl[[20,i,j]]=IMAGE[[20,i,j]]OTFl[[i,j]];]] 

For[i=l,  i<=32,  i++, 

For[j=l,  j<=32,  j++, 

IMAGEN 1  [[  1  ,i,j]]=IMAGE[[  1  ,i  ,j]]OTFn  1  [[i,j]] ; 
IMAGEN 1  [[2,i,j]]=IMAGE[[2,i,j]]OTFnl[[i,j]]; 
IMAGENl[[3,i,j]]=IMAGE[[3,i,j]]OTFnl[[i,j]]; 


IMAGENl[[20,i,j]]=IMAGE[[20,i,j]]OTFnl[[ij]];]] 

For[z=l,  z<=20,  z++, 

imagen  1  [[z]]=Abs[cinfft[IMAGEN  1  [[z]]]] ; 
imagep  1  [[z]]-Abs[cinfft[IMAGEPl  [[z]]]] ;] 

For[z=2,  z<=19,  z++, 

estimate[[z]]=image[[z]]  -  n  *  (imagenl[[z-l]]+imagepl[[z+l]]);] 

estimated  1  ]]=image[[  1  ]]  -  2n  *  imagepl[[l+l]]; 
estimate[[20]]=image[[20]]  -  2n  *  imagen  1  [[20-1]]; 

Figure  8.1:  Mathematica  implementation  of  the  Nearest  Neighbor  algorithm 
for  20  images. 


and  corresponding  restored  objects 


Figure  8.3:  Restoration  of  Mercury  spectral  radiance  curve  with  Nearest 
Neighbor  technique,  (a)  the  recorded  spectra,  and  the  reconstructed  spectra 
with  (b)  n=0.45. 

The  Nearest  Neighbor  technique  worked  for  the  above  application;  however,  it  is 
important  to  recognize  that  the  doublet  X  target  has  only  a  few  high  spatial  and  spectral 
frequencies.  The  computer  simulation  of  section  6.2  predicts  that  this  target  would  create 
little  blurring,  but  as  the  target  size  increases  and  the  emission  broadens  the  blurring 
becomes  stronger  and  the  reconstruction  more  challenging.  It  becomes  important  to  choose 
an  algorithm  which  takes  into  account  the  influence  of  the  whole  image  cube,  not  just  the 
adjacent  slices,  such  as  Inverse  Filtering. 
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8.3  Inverse  Filtering 

A  more  accurate  approach  to  remove  the  out-of-focus  information  requires  that  the 
full  contributions  of  all  of  the  observed  data  be  utilized.  The  most  direct  way  to  utilize  the 
entire  image  cube  is  to  perform  the  inverse  filter  operations  presented  in  section  4.2  and 

4.3.  0  =  ^1 

Starting  with  the  object,  image  and  point-spread  function  three-dimensional  data 
cubes,  recall  equation  4. 1 1 : 

i(x,y,z)  =  h(x,y,z)***  o(x,  y,  z) 

A  technique  used  in  representing  multi-dimensional  data  sets  is  lexicographical 
ordering.  Each  three-dimensional  matrix,  such  as  o(x,y,z),  can  be  converted  to  an  N3  x  1 
column  vector  by  using  a  stacking  operator  [Hunt  1977].  Applying  lexigraphical  ordering 
to  32x32x32  data  cubes  forms  vectors  with  32768  elements  each,  leaving  a  one 
dimensional  convolution: 

*'(v)  =  h(v)*  o{v )  (8.7) 

Taking  the  Fourier  transform  of  each  vector: 

I(y/)  =  H(Y)0(¥)  (8-8) 

The  object's  Fourier  spectrum  vector  0( y)  can  be  computer  from: 


CW  = 


m 

H{  yr) 


(8.9) 


An  inverse  Fourier  transform  is  applied  to  0(y/)  to  reconstruct  o(v),  which  can  be 
repartitioned  in  three  dimension. 

The  three-dimensional  matrix  calculation  has  been  reduced  to  a  scalar  division. 
Notice  that  consideration  must  be  given  to  the  possibility  of  dividing  by  zero.  Another 
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concern  is  that  the  regions  in  the  image  spectrum  where  the  OTF  values  are  very  small 
correspond  to  frequencies  that  are  mainly  noise.  Dividing  these  frequencies  by  the  low 
values  of  the  OTF  would  cause  an  over  amplification  of  noise.  Two  methods  are  generally 
used  to  deal  with  this:  apodization  and  applying  a  regularized  inverse  filter. 


8.3.1  Apodization  Inverse  Filter 

Apodization  is  a  technique  to  deal  with  dividing  the  image  spectrum  with  the  lowest 

values  of  the  OTF,  by  simply  setting  a  maximum  for  the  value  of  1  equal  to  a  constant 

H(W) 

a.  Values  greater  than  a  are  set  to  zero  [Erhardt  1985]. 


H(yr) 


1 


if  <  a 


H(ilf) 

0  if  >  a 


(8.10) 


An  apodized  inverse  filter  was  implemented  in  Mathematica.  The  code  is  presented 
in  Figure  8.4.  The  lexigraphical  ordering  function  is  performed  by  the  command  Flatten[ 
]  and  was  used  throughout  the  image  processing,  allowing  the  data  cubes  to  be 
"unwrapped".  The  function  Partition[  ]  is  used  after  the  processing  has  been  performed 
to  regenerate  the  three-dimensional  organization  of  the  data  cube.  Erhardt  used  this 
approach  with  a  theoretical  OTF.  However,  the  differences  between  experimental  and 
theoretical  OTFs  are  significant;  hence  better  results  were  obtained  using  the  experimental 
OTF  of  Figure  3.17.  Figure  8.5  contains  plots  of  lexicographical  vector  versions  of  the 


measured  3D  psf  and  OTF.  Figure  8.6  shows  the  results  obtained  from  applying  the 
apodized  inverse  filter  to  the  challenging  iris  target  and  its  computer  generated  16  pixel 
counterpart. 
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vector  element  #  v 

(a)  _ 


5000  10000  15000  20000  25000  30000  32768 

vector  element  #  v 

(b)  _ _ 


Figure  8.5:  List  Plot  of  the  flattened  (a)  psf.  (b)  OTF. 
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Figure  8.6:  Object  results  from  apodization  reconstruction  with  CG  images 
(left  column)  and  recorded  iris  images  (right  column). 
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8.3.2  Regularized  Inverse  Filter 

In  inverse  filtering  situations,  it  is  standard  practice  to  use  a  regularized  inverse 
filter  to  minimize  the  effects  of  noise  that  can  dominate  at  high  spatial  frequencies  [Sezan 
1992].  The  form  of  the  regularized  inverse  filter  for  the  scalar  equation  (8.9)  is: 


O(w)  =  - -  I  yr) 

l  HWH(yf)  +  a) 


(8.11) 


where  a  is  a  constant  whose  choice  is  based  on  the  signal-to-noise  ratio  of  the  data  and 

should  be  in  the  range  of  0.001  to  0.1.  In  the  experimental  reconstruction  a  value  of 

a=0.005  was  found  to  be  appropriate.  Notice  that  like  the  apodization  technique,  the  a  in 

the  denominator  of  equation  (8.11)  essentially  sets  a  threshold  maximum  of  _L_to  deal 

H(y/) 

with  the  low  values  of  the  OTF. 


A  regularized  inverse  filter  can  also  be  applied  in  the  three-dimensional  inverse 


equation: 


o(i,^o=H(^^cr'K^^o 


(8.12) 


0(£  yr,Q  =  [h(£  w£)*  H(£  Y’C)  +  «]"'  H(£  ^  0  (8-13) 


however  this  requires  taking  the  inverse  of  the  [H  H+a]  matrix,  a  task  which  can  be 
difficult  if  the  data  is  singular.  Further  discussion  can  be  found  below  in  section  8.3.3. 

The  scalar  regularized  inverse  filter  of  equation  (8.11)  was  implemented  in 
Mathematica  and  applied  it  to  the  32x32x32  recorded  image  cube  of  the  same  Mercury  X 
target.  The  code  is  printed  in  Figure  8.7. 
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flatimage=Flatten[image] ; 
flatpsf=Flatten[psf] ; 

FLATIMAGE=cfft[flatimage] ; 

FLATOTF=cfft[psf] ; 

a=.005; 

FLATOBJ=FLATIMAGE*Conjugate[FLATOTF]/((FLATOTF*Conjugate[FLATOTF])+a); 

flatobj=Abs[cinfft[FLATOB  J]] ; 

parobj=Partition[flatobj,  1024] ; 

3Dobj=Table[0.,{z,32},{x,32},{y,32}]; 

For[j=l.,  j<=32,  j++, 

3Dobj[[j]]=Partition[parobj[|j]],32];]; 

t3Dobj=Table[0.,{  z,32 } ,  { x,32 } ,  { y,32 }  ] ; 

For[j=l.,j<=32,j++, 

t3Dobj[[j]]=Transpose[RotateRight[Transpose[RotateRight[3Dobj[[j]],l  6]],  1 6]] ;] ; 

Figure  8.7:  Mathematica  code  for  the  regularized  inverse  filter. 

The  results  from  applying  the  regularized  inverse  filter  to  the  Mercury  X  target  are 

presented  in  Figure  8.8.  Pictured  in  Figure  8.8(a)  when  cx=0  the  regularized  inverse  filter 

defaults  to  standard  inverse  filter  of  equation  (8.9).  The  results  show  that  the  reconstruction 
is  ineffective  because  of  the  noise.  The  addition  of  the  regularized  inverse  filter  noise  term, 

a,  significantly  improved  the  image  quality.  The  image  sets  in  Figure  8.8  (b),  (c)  and  (d) 
show  the  reconstructed  object  sets  with  a  equal  to  0.1,  0.01  and  0.005  respectively.  A 

large  a  such  as  0.1  in  Figure  8.8  (b)  is  similar  to  a  low  pass  filter.  Setting  a  to  a  smaller 

number  increases  the  band  pass  for  better  reconstruction  but  also  increases  the  possibility 
of  noise.  Figure  8.9  shows  the  normalized  spectral  radiance  of  one  pixel  before  and  after 
processing. 
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(d)  g=0.005 


Figure  8.8:  Resulting  object  sets  after  regularized  inverse  reconstruction 
with  various  a. 


Continuing  the  exploration  of  the  spectral/spatial  relationship  this  algorithm  was 
applied  to  both  the  computer  generated  and  experimentally  measured  increasing  aperture 
doublet  targets  from  Figure  6.7.  The  results  are  pictured  in  Figures  8.10,  8.11  and  8.12. 


As  before,  each  2D  image  is  a  z (k)  vs.  y  cross-section  of  a  3D  data  cube  with  a 


corresponding  spectral  curve  plotted  below  it. 
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aperture:  50|im 


400jim  5 


'9®* 


(a)  measured  doublet  objects 


■L 


(b)  Regularized  inverse  reconstruction  with  q=0.005 
easured  and  regularized  inverse  reconstructed  Hg  doublet  images 


(a)  a=0.0001  too  much  noise 


(b)  a=0.001  noise  and  reconstruction  balanced  but  not  very  good 


(d)  a=0.01  low  noise  little  reconstruction _ 

Figure  8.12:  Results  of  regularized  inverse  filter  on  measured  400(im  and 
iris  images. 
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This  lexicographic  regularized  inverse  filter  worked  very  well  on  the  measured  X 
images  of  Figure  8.8  and  the  CG  data  of  Figure  8.7.  Even  the  long  32  pixel  aperture  of  the 
right  hand  column  is  reconstructed.  However,  it  was  ineffective  on  the  experimentally 
measured  image  cubes  of  the  increasing  apertures  shown  in  Figures  8.1 1  and  8.12.  Notice 
the  unpredicted  peaks  at  the  edges  of  the  spectral  curves  caused  by  both  the  apodization  in 
Figure  8.6,  and  regularized  inverse  filter  applied  to  measured  image  cubes.  These  peaks  are 
artifacts  of  an  invalid  assumption  during  lexigraphical  ordering  that  the  psf  is  shift  invariant 
and  circular.  The  limited  success  that  has  been  achieved  is  due  to  the  fact  that  the  targets 
have  been  spatially  and  spectrally  located  in  the  center  of  our  3D  data  cube,  similar  to  being 
padded  with  zeros.  However,  as  our  target  increased,  i.e.  the  iris  aperture,  the  circular 
approximation  is  seen  as  the  psf  wrapping  around  to  the  other  side  of  our  data  cube.  This 
problem  is  addressed  by  not  making  the  circular  assumption,  by  explicitly  performing  the  z 
dimension  matrix  multiplication  and  implementing  an  SVD  algorithm  which  also  addresses 
probable  difficulties  such  as  a  non-square  data  cube  (MxNxN)  and  singular  OTF  matrices. 

8.3.3  Singular  Value  Decomposition  Inverse  Filter 

DOIS  records  a  series  of  2D  spectral  images,  i , ,  i2,  — in.  These  are  saved  as  a  3D 
image  cube  image[[x,y,z]].  The  following  expands  equation  (8.2)  for  a  j= 4,  four  color 
target.  It  describes  each  spectral  image,  ij(x,y)  as  the  superposition  of  each  spectral  object, 
°k(x>y)  convolved  with  a  defocused  point  spread  function,  hjk(x,y)  [Mooney  1995], 

h(x>y)  =  *  *ox{x,y)  +  h,2(x,y)  *  *o2(x,y)  +  A,3(x,y)  *  *o3{x,y)  +  ^{x.y)  *  *o4{x,y ) 

i2(x,y)  =  h2l(x,y)**o1(x,y)  +  h22(x,y)*  *o2 (x,y)  +  h,3 (x, y )  *  *o3 (x, y)  +  h,4 (x, y)  *  *o4 (x, y) 
i3{x,y)  =  K{x,y)  *  *ox{x,y)  +  ^2(x,y)  *  *o2(x,y)  +  ^3(x,y)  *  *o3(x,y)  +  ^4(x,y)  *  *o4{x,y) 
i4(x,y)  =  h4]{x,y )  *  *ox{x,y)  +  h42{x,y )  *  *o2{x,y)  +  fc43(x,y)  *  *o3(x,y)  +  hu{xty)  *  *o4(x,y) 

(8.14) 

Following  the  theory  of  section  4.2  for  simulation  and  reconstruction,  2D  Fourier 
transforms  are  applied: 
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72(|,f)=HJI(|,f)O,(4.?)  +  HK(l.<)02(l.f)  +  ^(«.f)0,(^f)  +  ffJ4(4f)O4(«,f)(g  ]5) 
/!(4.f)  =  «J,(4.f)o,(4.f)+«K(4,f)o2(4,f)+HM(4,f)oJ(4.f)+»M(4.f)o4(4,‘f) 
/4(4,f)  =  H4,(4.m(4.?)+H«(4.f)o2(4.f)+H«(4.f)03(4.f)+ff44(4.f)o4(4.f) 

Rewritten  in  matrix  form: 


"A  (4.4) 

Mm  h12(|,4)  hb(4,4)  h„(4,4) 

1 

_o 

O 

_ 1 

h(U) 

h2](4,4)  hk(4.4)  he(4.4)  «»(4.4) 

o2(4.4) 

him 

h„(4.4)  h31(4,4)  Hj3(4.4)  h31(4,4) 

o3(4.4) 

him 

h„(4.4)  h<2(4,4)  h„(4.4)  ^(4,4) 

o.(4.4) 

This  equation  can  be  evaluated  for  each  spatial  frequency  (£,£).  Lets  take  a  closer 


look  at  the  H  matrix  : 


H  = 


H„(4,4)  Hl2(|,4)  Mm  ff,4(«.f) 
H21(|,4)  H22(4,4)  H^U)  H24(4.4) 
h31(4.4)  h32(|,4)  »3>(4,f)  H*(4.f) 
1 H«(10  Htl(U)  H„(U)  H„(U) 


(8.17) 


In  Chapter  4  great  lengths  were  taken  to  show  that  the  point  spread  function  is  shift 
invariant  (SIV).  Applied  here,  it  means  that  each  Hjk  can  be  expressed  by  the  Fourier 

C 

transform  of  a  point  spread  function  with  an  appropriate  amount  of  defocus,  HAz,  where 
Az  =  i  -  j,  i.e.  H0  is  infocus,  Hj  is  1  step  away  from  focus  and  H}  is  1  step  in  the  other 
direction. 


H  = 


tf„(«,4)  ff2(4,4)  H,(U) 

h,(4,4)  Ho (4.4)  H,(4. 4)  h2(4.4) 
h_2(4,4)  Mm  Ho  (4. 4)  H,(4,f) 
[h3(4, 4)  h2(4,4)  h,({,4)  h0(4,4)J 


(8.18) 


This  is  a  block  Toeplitz  matrix  formed  from  Fourier  transforming  images  of  point 
source  at  various  amounts  of  defocus.  Typically  when  considering  a  large  number  of 
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spectral  channels  it  will  become  a  sparse  band  diagonal  matrix  since  the  approaches 
zero  at  significant  amounts  of  defocus. 

The  task  is  to  invert  this  H  matrix,  however  it  can  have  various  problems  which 
cause  it  to  be  singular  and  difficult  to  invert  [Mooney  1995].  There  can  fail  to  be  a  unique 
solution  if  one  or  more  of  the  N  equations  is  a  linear  combination  of  the  others,  a  condition 
called  row  degeneracy.  Or  if  all  equations  contain  certain  variables  in  the  exact  same  linear 
combination,  inversion  can  fail  due  to  column  degeneracy.  A  set  of  equations  that  is 
degenerate  is  called  singular.  Singular  value  decomposition  is  a  technique  to  solve  for  a 
solution  to  the  inverse  matrix  with  singular  matrixes  and  in  fact  diagnose  the  singularities. 
This  is  the  method  of  choice  for  solving  most  linear  least  squares  problems  [Press  1988]. 

The  SVD  method  is  based  on  the  linear  algebra  theorem  that  any  MxN  matrix  H 
whose  number  of  rows  M  is  greater  than  or  equal  to  its  number  of  columns  N,  can  be 
written  as  the  product  of  an  MxN  column-orthogonal  matrix  U,  an  NxN  diagonal  matrix  W 
with  positive  or  zero  elements,  and  the  transpose  of  an  NxN  orthogonal  matrix  V. 

(8.19) 

The  matrices  U  and  V  are  each  orthogonal  in  the  sense  that  their  columns  are  orthonormal, 
where: 

T  V  =  1  (8.20) 
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Since  V  is  square,  it  is  also  row-orthonormal,  VVT=1.  The  inverse  of  the  diagonal  matrix 


W  is  easily  determined  by  replacing  the  diagonal  elements  with  1/w. 


Y’ll 

Wu 

_ 

1 

/»«. 

W„„ 

nn  j 

(8.21) 


The  decomposition  and  determination  of  U,  W,  and  VT  is  accomplished  by  a 
standard  routine  found  in  Numerical  Recipes.  In  Mathematica  the  function 
SingularValues[H]  outputs  three  matrices  U,  the  elements  of  W  and  V.  The  inverse  of 
H  is  found  from: 


H1 


=  V 


(8.22) 


The  elements  of  W  are  called  the  singular  values  of  the  system.  They  help  visualize 
the  missing  cone  of  our  data,  Figure  8.19.  The  noise  can  be  surpressed  by  adding, a 
regularization  filter  to  the  wu  values: 


1  _  w 
w  w2+al 


Substituting  this  SVD  into  the  inverse  reconstruction  equation: 

o(6f)  =  H(!,f)-,i(«) 


o(!,f)  =  \(U) 


diagl 


w 


m 


l(&0  + 


cr 


uT(£,C)  •!(£,£) 


(8.23) 


(8.24) 


This  SVD  Inverse  filtering  algorithm  was  implemented  in  Mathematica.  The  first 
step  was  to  form  the  OTF  matrix.  The  best  results  were  obtained  from  using  a  combined 
point  spread  function  which  is  an  average  of  the  experimentally  measured  and  theoretical 
psfs  found  in  Figure  8.13.  The  measured  psf  provides  information  about  the  influence  of 
aberrations  from  our  DOE,  while  the  theoretical  psf,  calculated  from  the  geometric  cone. 
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provides  the  very  small  values  which  fall  below  the  detectable  range  of  the  detector  and 
would  mistakenly  be  set  to  zero. 


object  reconstruction  (contrast  enhanced  to  show  detail). 

Looking  again  at  equation  (8.18)  one  might  notice  that  it  takes  seven  values  of 

Haz(^C)  to  form  the  four  z (k)  rows  in  the  H  matrix.  To  reconstruct  OTF  with  32  z(X) 

rows,  values  from  63  z(A.)  planes  are  required  as  shown  in  Figure  8.13.  This  provides 

images  of  a  point  source  at  focus  (plane  number  32)  and  31  defocus  steps  on  either  side  of 
focus.  Since  it  has  been  established  that  the  psf  is  shift  invariant  ,  a  "basis"  OTF  matrix  P 
is  formed  from  taking  the  2D  Fourier  transform  of  the  image  at  each  step  of  defocus  with 
the  procedure  listed  in  Figure  8.14. 


Figure  8.14:  Mathematica  code  to  calculate  the  "basis"  OTF,  P[[Az, £,£]]. 

The  OTF  matrix  is  formed  by  substituting  the  appropriate  defocus  Hata 
P[[Az,^]].  OTF  is  a  four-dimensional  matrix,  NxN  (32x32)  for  each  spatial  frequency 


This  code  takes  the  Fourier  Transform  of  the 
(32x,32y)  images  at  63  steps  of  defocus,  z, 
in  the  combined  psf  of  Figure  8.13,  then 
normalizes  the  results  to  form  the 

normalized  OTF,  P[[Az,£,Q] 


For [z=l, z  <  =  6  3  , z++, 

P[ [z]]=cfft [combinedpsf [ [z]  ]  ]  ; 
]  ; 

P  =  P /Max [Abs  [P]  ]  ; 


(£,,Q,  shown  in  Figure  8.15.  The  values  are  complex,  Re  +  ilm,  yet  the  standard  SVD 


code  is  written  for  real  elements.  However,  SVD  can  handle  non  square  matrices,  MxN,  so 
our  NxN  complex  matrix  is  converted  into  a  2NxN  real  matrix  where  the  first  N  rows  are 
the  real  values  and  the  remainder  are  the  imaginary.  It  is  this  matrix  which  will  be  inverted 
with  the  SVD  routine  and  then  rewritten  as  an  NxN  complex  inverse  matrix  using  the 
conversion  matrix  inim,  formed  in  Figure  8.16. 


OTF=Table[0.,{i,32>,{ j,32>, {  ^32>,{  £,32)],  "forms  a  4D  OTF  matrix” 

SessionTime [ ] 

For  [^=1/  %  <=32,  \  +  +  ,  "puts  the  proper  values  into  the  4D  OTF 

For[£  =  l,  C  <=32,  C  ++,  matrix  for  each  spatial  frequency  4,0" 

OTF [[££]]= 

{{Ptt324.cn,  Ptt314.cn,  Ptt304.cn,  ptt294.Cn,  ptt---4.CN,  piu.^]) 

{ P  [  [33 ,  C.O  ]  ,  P[[32,iCn.  P[[31,W1].  P[[30,wn,  Pft---4.CN,  P[[2,^N) 

{Ptt344.cn,  Ptt334.cn,  Ptt324.cn,  pff314.Cn,  p[[...4,Cn,  ptt34.0n 

{Ptt354.cn,  Ptt344.cn,  Ptt334.cn,  Ptt324.cn,  ?t t . . . 443 3 ,  Ptt44.cn} 

{Ptt364.cn,  Ptt354.cn,  Ptt344.cn,  Ptt334.cn,  Ptt...,Wii,  pt [544] 3 > 

{ p t [37 4, Cl ] ,  P[[364,cn,  P[[354,cn,  ptt344,C]],  Pit---4.CN,  Ptts-Wu 

{Pt [38  4.C] ] ,  P[[37  4.a],  Ptt364.cn,  P[[35  4,cn,  ptt---4.Cn,  Pin.Win 

{ p t [39 4, C] ) ,  Ptt384.cn,  Ptt37  4.cn,  Ptt3e4.cn,  Ptt---4.Cn,  Pt [8,441 3 } 

{p[[4o4,cn,  pt [39441 } ,  ptt384.CN,  Ptt37  4.C]],  Ptt-.-4.CN,  Ptt94,0n 

{Ptt4i4.cn,  p[[4o4,Cn,  Pt [39 4, Cl ] »  Ptt384.Cn,  Ptt---4.Cn,  p[t104>dn 

{P[ [42  4.C]  1 ,  P[[414,CN,  P[[404,CN,  Pt [394.C3 ] ,  Ptt---4.cn,  Pttn4.CN) 

{Ptt434.cn,  p[[424,C]),  pc [4i4.cn ,  Ptt404.cn,  Ptt---4.Cn,  Ptt124.Cn> 

{Pt [44  4.C3 ] ,  P[[434,CN,  P[[42  4,CN,  p[[4i4,C]],  Ptt.--4.Cn,  Ptti34.Cn> 

{pt [45 4, Ci ] ,  p[[444,cn,  Ptt434.cn,  Ptt424.cn,  ptt---4.Cn,  Ptti44.Cn> 

{Pt  [46  4,C>  ]  ,  P  [  [45  4,C>  > ,  P[[444,CN,  P[[434.CN,  Pft-.-4.CN,  P[  [154.0  ]> 

{P[[47  4,CN,  Pt  [46  4,C>  1  ,  Pt  [45  4,0  >  ,  P[[444,CN,  P[[---4.0>,  P[  [164.0  ]  > 

{Pt [484,0],  Pt [47  4.cn,  Pt [46  4.cn,  ptt454.Cn,  Ptt---4.Cn,  Ptti7  4.Cn> 

{p[ [49  4.cn,  Ptt484.cn,  Ptt474.cn,  Ptt4e4.cn,  ptt...4,Cn,  p[[i84,0n 


{Pt [6i4.cn,  Ptt604.cn,  Ptt594.cn,  Ptt584.cn,  pc t ■ 
{Pt [62  4.cn,  Pt [6i4.cn,  Ptt604.cn,  Ptt594.cn,  Ptt- 
{Pt [63  4.C1  ] ,  Ptt624.cn,  pt [614,0 ] ,  pi [604.cn,  Ptt- 

] ;  ]  ? 

SessionTime [ ] 

71630.88 

73560.42 


.4,Cn,  p[ [304,03) 

•  4,Cn,  Ptt3i4.cn> 

•  4,Cn,  p t [32 4.C3 ] > > 


Figure  8.15:  Mathematica  code  to  form  the  OTF  Matrix. 
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im=Table [0 . ,{x,32},{y,64}]; 
realim=IdentityMatrix [32] ; 

For [x=l , x<=32 , x++. 

For [y=l, y<=32 , y++ , 
im[ [ x , y ] ] =real im [ [x,y] ] ; 
im[ [x,y+32] ] =X*realim[ [x,y] ] ; ] 
;  ]  ; 


Li stDans ityPlot [Abs [ im] ] ; 


inim=PseudoXnverse [ im] ? 
ListDensityPlot [Abs [inim] ] ; 


The  SVD  algorithm  in  Numerical  Recipes 
and  Mathematica  only  works  on  real 
matrices  and  the  OTF  matrix  has  complex 
numbers.  However,  it  is  possible  to  rewrite 
the  square  NxN  complex  matrix  as  a  2NxN 
matrix,  called  ri,  with  the  real  values  in  the 
first  N  rows  and  the  imaginary  values  in  the 
lower  N  rows.  The  SVD  is  used  to 
decompose  and  invert  the  2NxN  ri  matrix. 

The  code  on  the  left  forms  a  2NxN  matrix 
inim  where  the  1st  N  rows  form  an  Identity 
matrix  and  the  bottom  rows  form  a  diagonal 
matrix  with  the  complex  value  i  along  the 
diagonal.  This  matrix  can  be  multiplied  by 
the  inverse  of  ri  to  reorder  the  inverse 
matrix  into  a  NxN  complex  inverse  matrix. 

H  I 


I 


Re[H] 


Im[H] 


SVD 


H 


-1 


S3  S3 

3  S' 

Pm  m 


Figure  8.16:  Mathematica  code  to  create  inim  which  converts  the  NxN 
Complex  matrix  into  an  2NxN  of  all  Real  values. 

The  code  in  Figure  8.17  creates  the  matrix  IMAGE  by  taking  the  Fourier  transform 
of  the  2D  spectral  images  at  each  z(k).  As  shown,  an  additional  dimension  was  added  to 

the  IMAGE  matrix,  m,  where  each  value  of  m  represents  a  different  target.  Since  the  OTF 
doesn't  change,  the  algorithm  can  be  applied  to  many  image  cubes  while  having  to  invert 


the  matrix  only  once.  In  fact,  once  a  proper  OTF  and  a  value  are  established,  the  inverse 
need  only  be  performed  once.  It  can  be  stored  in  memory  for  future  use. 
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“forms  a  matrix  with  Fourier  Transforms  of  6  measured  image  cubes  to  reconstruct  6  object 
cubes  simultaneously  ” 

I MAGE =T able [0. , {m, 6}, (z, 32} , {x, 32}, {y,32}] ; 

RE IMAGE = Table [0. , {m, 6} , {z,32}, {x,32}, {y,32}] ; 

For [m= 1 , m<=6,m++. 

For [z=l, z<=32, z++, 

IMAGE [[m,z]]=cfft [doublet image [ [m, z] ] ] ; 

For[x=l,x<=32,x++, 

For[y=l,y<=32,y++, 

REIMAGE [ [m, x, y , z ] ] = IMAGE [ [m , z , X , y } ] ; 3  ? 1 ; ] ; ]  ; _ 

Figure  8.17:  Mathematica  code  to  create  the  Fourier  space  matrix  IMAGE. 


INOTF=Table [0.,{x,32},{y,32}]; 
ri  =  Table [0., (x, 64 } , {y,32}]  ; 

s ingul a rvalue  s  = Table [0.,  { x , 3  2 } , {y,32}, { i , 3  2 } ] ; 

OBJ  =  Table [0 . , (m, 6 } , { z , 3 2 } , {  x ,  3  2  }  , {  y ,  3  2  }  ]  ; 
a=  .  0  9 

SessionTime [ ] 

For [x=l, x<=32 , x++ , 

For [y=l,y<=32,y++. 

For [k=l,k<=32,k++,  "form  the  2NxN  real  valued  matrix" 

ri [ [k] ] =Re [ OTF [ [x,y,k] ] ] ; 

ri [ [k+32] ] =Im[OTF [ [x, y , k] ] ] ; ] i 

[u,md,v}=SingularValues  [ri]  ;  "use  the  SVD  command  to  find  the 
singularvalues  [  [x,y]  ]  =md;  three  matricies  U,W  &  VT" 

w=md / (md*2+a*2) ; 

inri=Transpose [v] . DiagonalMatrix [w] .u;  "find  the  inverse" 

INOTF  =  inri  .  inim;  "reorder  into  the  complex  inverse  matrix" 

For [m=l,m<=6,m++,  "Find  the  object  from  O^'1!" 

OBJ [ [m, x,y] ] =INOTF . REIMAGE [ [m, x,y] ] ; ] ; 

1  i 

Print [x] ; 

1  ; 

SessionTime [ ] 

71630.88 

73563.42 

{Max[singularvalues],Min[singularvalues]} 

(14.5595,  0.0000512428}  _ _ _ 

Figure  8.18:  Mathematica  code  finds  the  SVD  matrices,  inverts  the  OTF  and 
finds  the  reconstructed  OBJ  matrix. 

The  singular  value  decomposition,  inverse  and  object  reconstruction  are  performed 

with  the  code  listed  in  Figure  8.18.  The  matrix  ri  is  the  2NxN  matrix  of  real  values.  The 

function  SingularValues[ri]  outputs  the  three  matrices  u,  md,  and  v.  The  singular  values 
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md  are  stored  for  future  investigation  and  the  a  noise  suppression  is  applied.  The  inverse 

is  stored  in  INOTF  then  multiplied  by  the  IMAGE  resulting  in  OBJ.  This  is  then  inverse 
Fourier  transformed  to  reconstruct  the  object  cube  in  the  program  listed  in  Figure  8.19. 

"Reorder  the  matrices"  “ 

REOBJ=Table [0 . ,  (n,  6 }  ,  { i f 3 2 } , { j , 3 2 > , {k, 3 2 } ] ; 

ob j ect  =  Table [  0  .  ,  {n,  6 } , { i , 3 2 > , { j , 3 2 > , {k, 3  2 }  ]  ; 

For [z=l, z  <  =  3  2  , z  +  + , 

For[x=l,x<=32,x++, 

For [y=l,y<=32 , y++. 

For [m=l,m<=6,m++, 

REOBJt [m, z , x, y] ] =OBJ[ [m,x,y, z] 

"Take  the  inverse  Fourier  Transform  of  OBJ  to  find  the  object  cubes" 

For [m=l,m<=6,m++. 

For [z  =  l, z  <  =  3  2  # z  +  +  , 

object [ [m,  z  ] ] =Abs [cinf f t [REOBJ[ [m, z] ] ] ] ; 

] ; Print [m] ] ; _ 

Figure  8.19.  Muthcmutico,  code  takes  the  inverse  Fourier  transform  of  the 
OBJ  matrix. 

The  singular  values,  of  the  OTF  matrix  are  plotted  in  Figure  8.20.  Each 

image  is  a  density  plot  of  the  spatial  frequencies  contained  at  a  given  spectral  frequency. 
Notice  that  there  is  a  region  of  black  pixels  in  the  center  of  each  image  which  grows  with 
increasing  spectral  frequencies.  This  represents  low  frequency  spatial/spectral  information 
that  isn't  imaged  by  DOIS  or  recoverable  with  processing.  Plotted  in  Figure  8.21  as  a 
cross-section  of  each  spectral  frequency,  the  frequently  discussed  missing  cone  which 
plagues  such  imaging  systems  is  seen  [Barrett  1981,  Mooney  1995,  Descour  1995], 

Singular  value  decomposition  inverse  reconstruction  provides  a  unique  and  valuable 
opportunity  to  quantify  the  missing  cone.  The  missing  spatial  and  spectral  frequencies  can 
be  determined  from  the  design  parameters  to  predict  the  cut  off  frequencies  and  used  as  a 
figure  of  merit  when  designing  and  comparing  future  imaging  spectrometers. 
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(a)  wnn(£,0  increasing  spectral  frequencies — : 


(b)  increased  contrast  on  1st  10  wnn(£,Q  increasing  spectral  frequencies- 


(c)  the  10  images  of  (b)  only  larger 


Figure  8.20:  Singular  Values  w„(£,0  Matrix  plots. 


Figure  8.21:  Cross-section  shows  missing  cone. 


8. 3.3.1  SVD  results 


The  following  figures  present  results  from  applying  the  SVD  object  reconstruction 
algorithm  from  equation  8.24  to  computer  generated  (CG)  and  experimentally  measured 

images  cubes,  (32x,  32y,  32z(X))  in  size.  Figure  8.23  depicts  the  results  from 

reconstructing  the  CG  object  cubes  of  increasing  square  apertures  in  front  of  a  Hg  doublet 
that  were  generated  in  section  6.2.2.  As  expected,  the  doublet  is  indeed  reconstructed  using 
the  SVD.  Looking  closer  at  the  results,  Figures  8.24  and  8.25  picture  the  same  CG  data 
sets  but  displayed  as  series  of  spectral  images.  The  first  is  from  the  8x8  square  aperture, 
and  the  second  from  the  16x16  square.  Objects  were  reconstructed  with  the  regularization 

filter  value  a  set  to  0.5,  0.1  and  0.05  from  equation  8.23.  Notice  that  only  30  images  are 
printed  for  each  image  series,  for  convenience  the  1st  and  32nd  were  omitted  from  the 
figure.  Figure  8.22  lists  the  location  values  z(X)  corresponding  to  the  (32x,  32y)  images 
pictured  in  the  series.  Recall  that  each  step  represents  dz~0. 125mm  and  spectral  step  size 
dX~  25nm  in  both  the  CG  and  experimental  cubes.  The  entire  32  step  series  records  an 
8nm  band,  a  resolution  that  rivals  any  spectral  filter. 
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Figure  8.22:  Table  of  the  z  values  for  the  following  output  series. 


In  (a)  of  Figure  8.24  and  8.25,  the  square  aperture  emits  at  locations  14  and  20. 
The  effect  of  DOIS  imaging  is  seen  in  (b):  the  edges  of  the  squares  are  blurred,  and  it 
inaccurately  measures  the  presence  of  the  square  in  more  than  the  two  locations,  the  image 

is  now  blurred  over  most  z  locations.  The  SVD  Inverse  filter  algorithm  was  applied  with  a 
equal  to  0.5,  0.1  and  0.05  with  the  results  displayed  in  (c),  (d)  and  (e)  respectively.  The 
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reconstruction  does  reduce  the  occurrence  of  the  square  in  the  inaccurate  locations,  and  the 
square  is  sharpened.  As  a  is  made  lower,  the  reconstruction  is  better;  however,  noise  and 

ghost  images  are  introduced  which  aren't  present  in  the  (b)  image  cube. 

Interestingly,  the  reconstructed  object  squares  appear  to  be  shifted  by  two  z 

locations,  dA.^0.5  nm.  The  cause  of  this  shift  is  indeterminate,  an  investigation  into  this  is 

matter  is  recommended  as  a  suggestion  for  further  work. 

The  SVD  algorithm  was  applied  to  four  different  experimentally  measured  image 
cubes.  Each  target  was  assembled  to  demonstrate  a  variety  of  spectral/spatial 
characteristics.  For  consistency  Figure  8.26  shows  that  the  SVD  reconstmction  works  well 
on  the  experimentally  measured  577/579  nm  Mercury  doublet  with  the  X  aperture.  The 

noise  regularization  value  a  was  set  to  0.5,  0.1  and  0.05  for  (b),  (c)  and  (d)  respectively  in 
the  figures  of  all  four  targets. 

The  second  target  combines  a  542nm  GreNe  laser  point  source  imaged  along  side 
the  546  nm  Hg  lamp  with  a  trapezoidal  aperture.  The  recorded  image  cube  and 
reconstructed  object  cubes  are  shown  in  Figure  8.27.  The  SVD  reconstruction  improved 
both  the  spectral  and  spatial  quality  of  the  images  without  the  artifacts  demonstrated  with 
the  other  algorithms. 

To  touch  on  larger  targets,  an  almost  full  FOV  rectangle  aperture  was  placed  in 
front  of  the  Hg  doublet  in  Figure  8.28.  Although  the  2nm  doublet  can't  be  resolved,  DOIS 
measures  it  to  be  a  line  with  about  4  nm  FWHM  using  the  SVD  reconstruction,  and  the 
spatial  information  of  the  rectangle  aperture  is  reconstructed.  Finally,  the  mercury 
illuminated  iris  image  cube  was  processed  with  the  SVD.  The  results  are  shown  in  Figure 
8.29,  notice  that  a  "hot  spot"  on  the  bulb  has  been  reconstructed. 
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(c)  SVD  reconstruction  with  a=.5 
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(d)  SVD  reconstruction  with  oc=.l 
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J  _ (e)  SVD  reconstruction  with  a-.Q5 _ _ 

Figure  8.24:  SVD  reconstructed  object  series  from  the  CG  8x8  pixel 
doublet  aperture,  previously  depicted  as  cross-sections  in  8.23. 
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(a)  CG  doublet  square(16x!6pixels)  target  at  #13&  20 


(b)  CG  image  set  with  measured  psf 


(c)  SVD  reconstruction  with  oc=0.5 


(d)  SVD  reconstruction  with  a=.l 
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1 _  (e)  SVD  reconstruction  with  a=.Q5 _ 

Figure  8.25:  SVD  reconstructed  object  series  from  the  CG  16  xl6  pixel 
doublet  aperture,  previously  depicted  as  cross-sections  in  8.23. 
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Figure  8.26:  SVD  reconstructed  object  series  from  the  experimentally 
measured  Hg  577/579  nm  X. 
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(a)  experimentally  measured  images 
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| _ (d)  SYD  reconstruction  q=75 _ 

Figure  8.27:  SVD  reconstructed  object  series  from  the  experimentally 
measured  542  nm  GreNe  point-source  with  a  trapezoidal  546  nm  Mercury 
target. 
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Figure  8.28:  SVD  reconstructed  object  series  from  the  experimentally 
measured  Mercury  Rectangle. 


(notice  poor  reconstruction  when  the  target  extends  to  the  edges) 
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(b)  SVD  reconstruction  a=1.5 
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rigure  o.^y:  o\u  reconstructed  ooject  series  irom  the  experiments 
measured  Iris  illuminated  by  a  Mercury  bulb  (notice  the  hot  spot). 
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8.4  Constrained  Iterative  Deconvolution 

The  SVD  Inverse  Reconstruction  algorithm  with  a  regularized  filter  for  noise 
suppression  provides  a  linear  restoration.  One  drawback  of  such  linear  methods  is  that 
negative  ripples  build  up  around  strong  features  from  the  low  values  of  the  OTF  and  no 
reconstruction  is  possible  beyond  the  OTF  cutoff.  Nonlinear  restoration  methods  that 
constrain  the  result  to  be  positive  can  eliminate  this  problem  while  using  positivity  and 
other  a  priori  information  to  provide  a  better  restoration.  The  drawback  is  that  these 
methods  are  invariably  iterative  and  hence  quite  computationally  expensive  in  three 
dimensions.  An  additional  challenge  is  to  develop  a  relaxation  parameter  optimized  for  a 
particular  application  and  image  cube. 

The  strategy  is  to  develop  a  positively  constrained  solution  or  guess  object  cube  go, 
that  when  convolved  with  the  known  smearing  function  h  will  regenerate  the  observed 
image  cube  i..  The  pixel-by-pixel  differences  between  the  convolved  guess  and  the 
observed  data  are  used  to  update  the  guess. 

gik  -  gok  *  h 

gok  1  =go  +y(i-gi  )  (8.25) 

if  gok+x  <0  then  gok+l  =  0 
k  =  k  +  l 

This  is  the  Jansson-vanCittert  constrained  iterative  deconvolution  algorithm  with  the 
relaxation  parameter  [Jansson  1984]: 

r  =  r„[l-2|go*-!|]  (8.26) 

which  applies  the  a  priori  information  that  the  intensity  is  a  positive  value  and  that  the 
spectra  is  a  transmission  function  with  a  maximum  value  of  1.  The  best  results  were 
obtained  with  an  initial  guess  equal  to  the  regularized  inverse  filtered  data  set  and  a 

weighted  relaxation  parameter,  y. 
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The  vanCittert  iterative  technique  with  a  Jansson  relaxation  parameter  of  equation 
(8.25)  was  implemented  in  Mathematica  and  is  listed  in  Figure  8.30.  The  algorithm  was 
applied  it  to  the  X  Mercury  doublet.  The  results  after  10  and  20  iterations  are  pictured  in 
Figure  8.31.  The  spectral  radiance  curve  plotted  in  Figure  8.32  shows  that  the  doublet  is 
better  resolved  with  every  iteration. 


gama=flatmimage/Max[flatmimage] ; 
gobj=regularizedobj ; 

Print]  "=start  time"SessionTime[]];  Print["=original  maxobj"Max[gobj]]; 

Do[ 

GOB  J=cflatfft[gobj] ; 
gi=Abs[cflatinfft[GOBJ*FLATOTF]] ; 
diff=flatmimage-gi ; 
gobj=gobj+gama*diff; 
g°bj=(gobj + Abs  [gobj]  )/2 ; 
objnum[[c]]-gobj; 

Print["iteration"];Print[c];  Print["=maxobj"Max[gobj]];  Print["=maxdiff'Max[diff]]; 
Print["=mindiff'Min[diff]l ;  Print["time"SessionTime[]] ; 

_ ,{c,20}] _ _ _ 

Figure  8.30:  Mathematica  code  for  the  Jansson-vanCittert  algorithm. 
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Figure  8.31:  Resulting  images  after  Jansson-vanCittert  reconstruction  with 
various  iterations,  (a)  original  recorded  data,  (b)  n=10.  (c)  n=20. 
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The  Jansson-vanCittert  algorithm  didn't  work  for  large  targets  or  data  with 
excessive  noise.  Various  relaxation  functions  were  implemented  and  several  hundred 
iterations  ran  without  ever  reaching  convergence.  This  inspired  the  addition  of  another  step 
to  the  algorithm.  The  difference  term  is  inverse  filtered  with  a  regularized  inverse  filter  as 
described  in  equation  (8.27).  This  achieved  much  better  results  that  converge  within  a  few 
iterations. 


gik  =gok*h 
diff  =  i  ~  gik 
DIFF  =  Fourier[diff] 


regulardiff  =  Abs 


InverseFourien 


( 


DIFF] 


H* 


H*H  +  a 


(8.27) 


gok+i  =  gok  +  regulardiff 
if  gok+l  <  0  then  gok+l  =  0 
k  =  k  + 1 

The  regularized  inverse  filtered  vanCittert  algorithm  of  equation  (8.27)  was 
implemented  in  Mathematica  with  the  code  listed  in  Figure  8.33.  It  was  applied  to  CG 
doublet  and  experimentally  measured  iris  images.  The  rms  difference  between  the  recorded 
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images  and  the  guess  images  is  plotted  against  the  number  of  iteration  in  Figure  8.34;  they 
show  convergence  in  less  than  5  iterations. 

start  with  measured  image 
image=flatmimage/Max[flatmimage] ; 
gobpimage; 
a=.01 

Print["=start  time"SessionTime[]]; 

Do[  gimage=Abs[cflatinfft[cflatfft[gobj]  *FLATOTF]] ; 
gimage=gimage/Max  [gimage] ; 
diff=image-gimage; 
diff=diff-Min  [diff] ; 
n=Ceiling[(c)]; 

REGULARDIFF=cflatfft[diff]*Conjugate[FLATOTF]/((FLATOTF*Conjugate[FLATOTF])+a); 
regulardifff  [n]]= Abs  [cflatinfft[REGUL  ARDIFF]  ] ; 
gobj=gobj+regulardiff[[n]]; 
rms[[c]]=Mean[Abs[regulardiff[[n]]]]; 
objnum[[n]]=gobj; 

Print["iteration"c];  Print["=rms"rms[[c]]]; 

,{c,20}] 

Print["time"SessionTime[]]; 

For[i=l.,  i<=l,  i++, 
par[[i]]=Partibon[reblurredobject,  1 024] ; 

For[j=l.,  j<=32,  j++, 
obj[[i,j]]=Partition[par[[i)j]],32]; 

obj  [[i  j]]=Transpose[RotateRight[Transpose[RotateRight[obj[[i,j]] ,  1 6]],  1 6]] ;] ;] ; 

For[m=l,  m<=l,  m++, 

For[z=l,  z<=32,  z++, 
trans[[z]]=Transpose[obj[[m,z]]]; 
h[  [m,z]]=trans  [[z,  16]]; 
w[[m,z]]=obj  [[m,z,  16]];]; 
h[[m]]=Transpose[h[[m]]] ; 

w[[m]]=Transpose[w[[m]]];]; _ 

Figure  8.33:  Mathematica  code  for  regularized  inverse  filtered  vanCittert 
iterative  reconstruction. 


The  results  are  depicted  in  Figures  8.35,  8.36  and  8.37.  The  objects'  spectral  and 
spatial  information  is  reconstructed,  including  the  lamp  hot  spot  and  spectral  doublet  in  the 
iris  data.  Although  this  algorithm  does  look  promising,  it  requires  the  regularized  inverse 
filter  which  was  shown  to  make  approximations  that  aren't  always  appropriate.  See  the 
wrap  around  problems  in  the  iris  plots  which  result  in  improper  peaks  in  the  tails  of  the 
spectral  plots.  It  was  also  found  that  the  algorithm  had  to  manipulated  to  fit  the  target.  This 
aside,  the  constrained  iterative  technique  could  still  be  a  solution  to  the  deconvolution 
problem  for  certain  applications,  particularly  if  DOIS  were  required  to  image  the  same  or 
similar  targets  repeatedly. 
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Figure  8.37:  Reconstructed  iris  objects  after  10  &  20  iterations. 
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9.  CONCLUSION 


The  Diffractive  Optic  Image  Spectrometer  (DOIS)  is  a  viable  multi-spectral  imager. 

This  dissertation  demonstrated  and  proved  the  concept  of  spectral  imaging  with  a 
diffractive  optic.  More  importantly,  the  theory  of  spectral  imaging  with  a  diffractive  optic, 
referred  to  as  diffractive  spectral  sectioning,  was  developed  and  an  appropriate  digital 
image  restoration  algorithm  was  derived  to  deconvolve  the  DOE  transfer  function  from  the 
recorded  image  cube,  providing  a  high  resolution  spectral/spatial  object  cube. 

The  spectral  resolving  power  of  288  (KjAk  =  577nm/2  nm)  demonstrated  by  the 

DOIS  prototype  can  be  utilized  to  predict  capability  of  future  diffractive  optic  image 
spectrometer  designs  applied  in  any  spectral  range;  UV,  visible  or  IR. 

I  conclude  that  the  best  image  spectrometer  designed  with  a  diffractive  optic  would 
include  a  kinoform;  a  high  efficiency,  multi-level  DOE  with  a  refractive  surface  to  balance 
the  spherical  and  off  axis  aberrations.  An  intermediate  zoom  lens  should  be  used  to  perform 
the  spectral  scanning  while  maintaining  a  constant  magnification.  The  most  accurate  and 
powerful  reconstruction  algorithm  is  the  Singular  Value  Decomposition  Inverse  Filter  with 
a  regularization  term  for  noise  suppression. 

9.1  DOIS  Advantages 

DOIS  has  several  advantages  over  current  image  spectrometers.  It  incorporates  a 
simple  one  axis  translation  on  a  rugged  platform  making  it  insensitive  to  vibrations  which 
limit  Fourier  Transform  spectrometers.  DOIS  is  programmable,  providing  single  spectra, 
narrowband  or  full  spectrum  image  cubes.  It  can  provide  coarse  or  fine  spectral  resolution 
by  choosing  the  stepping  increment  and  an  object  reconstruction  algorithm  at  various  levels 
of  computational  expense. 
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DOIS  uses  mainly  off  the  shelf  components.  The  DOE  fabrication  does  have  a  large 
first  time  expense  to  generate  a  master;  however,  multiple  copies  can  be  replicated  at 
minimal  expense.  The  design  is  not  limited  by  the  availability  of  materials  like  conventional 
thin  film  spectral  filters.  Once  built  for  one  application,  additional  wavelengths  within  a 
wide  spectral  range  can  be  viewed  for  other  applications  without  changing  components. 

DOIS  solves  a  common  problem  associated  with  spectral  filters;  the  central 
wavelength  of  a  filter's  bandpass  can  shift  due  to  environmental  factors  such  as 
temperature.  This  can  be  corrected  for  in  a  DOIS  system  with  a  simple  change  in  position, 
making  on-board  calibration  and  realignment  possible. 

DOIS  provides  enough  spectral  and  spatial  image  quality  without  post-detection 
processing  that  there  are  applications  where  the  recorded  image  cube  can  be  utilized  to 
represent  the  target.  Contrary  to  computer-tomography  approaches,  even  when  object 
reconstruction  processing  is  required  the  pre-processed  images  are  at  least  recognizable  so 
that  the  operator  can  have  confidence  while  recording  the  data. 

The  DOIS  spectrometry  feature  can  be  added  to  existing  camera  systems  with  a 
simple  lens/mount  replacement,  providing  additional  information  for  the  difficult  task  of 
target  identification.  DOIS  can  also  be  a  cost  effective  solution  to  spectroscopy  applications 
where  imaging  isn’t  required  but  would  serve  as  a  bonus.  The  imaging  will  help  in 
minimizing  misalignment  and  improve  tracking  of  moving  targets.  The  additional  pixels 
within  the  field  of  view  could  act  as  simultaneous  data  measurements  that  can  be  averaged 
to  improve  accuracy  and  reduce  noise. 

As  with  most  imaging  spectrometers  the  scanning  mechanism  will  limit  use  in 
applications  with  short  lifetimes,  however  as  presented  below  in  section  9.2.3,  a  DOIS 
system  can  be  designed  with  no  moving  parts,  recording  the  entire  image  cube  in  a  single 
"snapshot". 
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9.2  Suggestions  for  future  work 


During  this  dissertation  several  ideas  were  conceived  for  continued  research  with 
diffractive  optic  imaging  spectrometers  and  are  presented  below.  Additionally,  it  was 
discussed  in  section  6.2.2  that  a  field  stop  should  be  added  to  the  optical  train  to  aid  in 
reconstructing  targets  on  the  edge  of  the  field  of  view  and  that  an  investigation  should  be 
performed  into  the  cause  of  a  spectral  shift  in  the  reconstmcted  data  in  Figure  6.7. 

9.2.1  Dual  Waveband  design  using  multiple  orders 

While  working  with  the  2  level  zone  plate,  5  separate  diffractive  orders  were 

observed  by  scanning  the  detector  closer  to  the  DOE.  Each  order  is  located  at  — ,  where  f 

m 

is  the  focal  length  of  the  first  diffractive  order  and  m  is  the  diffractive  order  number.  While 
viewing  wideband  emitting  targets,  the  orders  started  to  superimpose  one  another.  Why  not 
use  this  effect  to  design  a  multi-band  spectral  imager,  where  the  DOE  images  far  IR 
radiation  in  the  first  order  and  the  mid  IR  in  the  second  order? 


Figure  9.1:  Two  waveband  design,  mid  IR  and  far  IR,  with  X0=7.3(am 
spectral  r\  for  each  order,  and  the  two  wavelengths  at  each  focal  position. 

A  16  level  DOE  designed  at  A,d=7.3pm  images  8  to  12pm  far  IR  radiation  in  its  first 

order  and  3  to  5pm  mid  IR  in  its  second  diffractive  order.  Figure  9.1  shows  the  1st  and  2nd 
order  spectral  diffraction  efficiency  for  a  Xd=7.3pm  DOE.  The  two  wavelengths  that  are 
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imaged  at  each  focal  position  are  listed  in  Figure  9.1(b).  One  can  think  of  various  detector 
configurations,  as  well  as  an  array  of  both  mid  and  far  IR  detectors. 

This  multi-order  concept  was  demonstrated  with  the  DOIS  visible  prototype.  The 
second  order  image  of  a  633  nm  HeNe  laser  came  to  focus  at  the  same  plane  as  the  third 
order  image  of  a  422  nm  monochrometer  slit. 

9.2.2  Programmable  DOE  with  an  SLM,  variable  focus  device 

Depicted  in  Figure  9.2,  an  electronically  controlled  variable  focus  diffractive  lens  or 
spatial  light  modulator  (SLM),  can  replace  the  DOE  in  the  diffractive  optic  image 
spectrometer  to  create  a  DOIS  with  no  moving  parts.  The  components  are  mounted  with  a 
fixed  separation.  The  SLM  will  select  the  image  spectra  by  electronically  changing  the  DOE 
ring  spacing,  which  will  change  the  focal  length  [Domash  1996]. 


Figure  9.2:  Design  for  a  static  DOIS  with  a  controllable,  variable  focus 
lens. 

9.2.3  Array  of  DOEs 

Finally,  an  image  spectrometer  can  be  designed  with  diffractive  optics  and  no 
moving  parts  to  record  the  entire  3D  image  cube  simultaneously.  Figure  9.3  is  the 
schematic  of  a  DOIS  "snapshot"  technique  which  records  a  set  of  images  at  predetermined 
spectral  lines.  Each  lens  of  the  array  is  fabricated  to  have  the  same  design  focal  length  but  at 
different  wavelengths.  The  detector  array  is  segmented  into  a  region  for  each  wavelength  or 


an  array  of  detector  arrays  could  be  assembled.  This  configuration  eliminates  scanning  and 
magnification  differences.  However  unlike  a  segmented  detector  with  an  array  of  filters,  if 
the  temperature  shifts  or  requirements  change  the  array  can  easily  be  repositioned  at  a 
different  distance  to  detect  another  set  of  spectra.  Additionally,  by  fabricating  one  lens  to  be 
refractive,  such  as  a  Fresnel  Lens,  a  conventional  full-band  image  of  the  entire  scene  can  be 
recorded  as  well. 


Figure  9.3:  Design  for  a  static  DOIS  with  a  DOE  lens  array. 


In  conclusion,  DOIS,  the  diffractive  optic  imaging  spectrometer  presented  in  this 
dissertation  is  a  practical,  high  resolution,  compact,  economical,  rugged,  programmable, 
multi-spectral  imager.  With  proper  selection  of  the  DOE  substrate  and  detector  material,  it 
can  be  built  to  operate  at  ultraviolet,  visible  or  infrared  wavelengths  for  applications  in 
surveillance,  remote  sensing,  law  enforcement,  environmental  monitoring,  laser 
communications,  medical  imaging  and  laser  counter  intelligence. 
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